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Abstract  of  Thesis  Presented  to  the  Graduate  School 
of  the  University  of  Florida  in  Partial  Fulfillment  of  the 
Requirements  for  the  Degree  of  Doctor  of  Philosophy 

DRAG  MEASUREMENT  ON  A SUSPENDED  SPHERE  AND  ITS  APPLICATION 

TO  CORROSIVE  GAS  VISCOSITY  MEASUREMENT 

By 

Ratan  Kxomar 
December  1992 

Chairman:  Dr.  Samim  Anghaie 

Major  Department:  Nuclear  Engineering  Sciences 

The  estimation  of  drag  force  exerted  by  fluid  flow  over 
a sphere  is  a fundamental  problem  and  has  received  attention 
both  for  bounded  and  unbounded  flow.  An  offshoot  of  this 
study  on  a sedimenting  sphere  in  a cylindrical  tube  has 
provided  an  effective  means  to  ascertain  viscosity  of 
fluids.  However,  low-drag  measurement  difficulties  and  other 
practical  impediments  had  restricted  the  sedimentation  study 
to  the  determination  of  liquid  viscosity  only.  The  current 
work  incorporates  an  accurate  microbalance  to  obtain  the 
drag  force  exerted  by  gas  flow  over  a suspended  sphere.  The 
drag  measurements  are  carried  out  by  an  ultra  precise 
thermal  microbalance  having  a sensitivity  of  0.1  microgram. 
The  results  obtained  have  been  extended  for  viscosity 
measurements.  The  experimental  findings  provide  enough 
impetus  to  propose  a methodology  that  can  assist  in  the 
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development  of  a drag  based  viscometer  with  good  potential 
to  handle  corrosive  gases. 

In  conjunction  with  the  experimentation,  a theoretical 
investigation  afforded  through  dimensional  analysis  and 
nvimerical  studies  has  been  carried  out.'^The  examined 
phenomenon  is  the  flow  of  a fully-developed  gas  in  a 
cylinder  over  a centrally  suspended  sphere.  Previous  studies 
in  this  area  had  come  up  with  correlations  that  established 
association  between  the  coefficient  of  drag  on  a sphere  in 
unbounded  flow  and  cylindrically  confined  flow.  These 
relationships,  however,  do  not  corroborate  the  experimental 
data  well.  This  prompts  the  developement  of  a new 
correlation  that  can  be  explained  physically. 

Finally  an  attempt  has  been  made  to  explain  the 
observed  experimental  results  by  computational  methods. 
Although  the  flow  geometry  is  fundamental,  it  is  difficult 
to  numerically  test  it  for  drag  measurements.  Two  different 
approaches  have  been  taken  in  the  computational  fluid 
dynamics  study.  A body-fitted  coordinate  system  using  the 
streamfunction-vorticity  method  and  a primitive  variable 
technique  in  cylindrical  coordinates  were  utilized.  The 
results  show  that  the  streamfunction-vorticity  method  runs 
into  some  difficulty  with  the  body-fitted  coordinate  system. 
The  primitive  variable  approach  appears  more  robust  and  can 
be  recommended  for  future  work  with  the  body-fitted 
coordinate . 
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CHAPTER  1 
INTRODUCTION 

1 . 1 Motivation 

Fluid  flow  in  nuclear  systems  provides  complex 
situations  which  are  difficult  to  handle  experimentally  and 
theoretically.  Several  studies  on  advanced  nuclear  reactors 
for  space  application  have  been  done  but  their 
implementation  on  a practical  scale  is  still  a serious 
challenge.  One  such  proposed  nuclear  system  for  space  power 
generation  is  shown  in  Figure  1-1.  Uranium  fuel  in  droplet 
form  is  sprayed  into  the  reactor  which  generates  about  750 
megawatts  of  thermal  energy.  In  order  to  simulate  and  test 
this  idea  from  a thermalhydraulic  point  of  view, 
thermophysical  properties,  like  viscosity,  of  the  working 
fluid  and  correlations  for  fluid  flow  around  spherical 
bodies  need  to  be  ascertained.  Highly  reacting  gases  can 
have  certain  unresearched  effects  on  the  measuring  devices. 
As  a result,  viscometers  that  are  exposed  to  corrosive  gases 
may  provide  unrealistic  data.  This  motivates  the  current 
study  to  propose  a viscometer  that  is  not  exposed  to  the 
unknown  environment  and  has  the  attractiveness  of  a precise 
and  sensitive  measuring  device.  Another  interesting  feature 
of  fluid  flow  around  spherical  bodies  enclosed  by 
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Figure  1.1  Schematic  of  200  MWe  uranium  droplet  fuel,  superheated 

lithium  reactor  with  mhd  generator 
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walls  is  also  discussed  in  this  work.  Past  studies  on  pebble 
bed  nuclear  reactors  stimulated  interest  in  turbulent  flow 
over  cylindrically  surrounded  spheres.  Theoretical  work,  in 
this  area,  has  been  performed  for  very  low  flow  rates.  In 
the  intermediate  laminar  region,  reliable  data  have  been 
generated  by  interpolating  between  the  turbulent  and  low 
flow  rate  regimes.  The  present  study  endeavors  to  highlight 
some  aspects  of  laminar  fluid  flow  over  a rigid  sphere  in  a 
tube  by  experimentation  and  analysis. 

1.2  Discussion  on  Fluid  Drag  Force  and  Viscosity 
When  an  object  is  placed  within  a flowing  fluid  stream, 
the  latter  experiences  a change  in  its  momentum  and  a force 
is  exerted  on  the  object.  This  force  owes  its  existence  to 
the  viscosity  of  the  fluid  and  to  the  pressure  distribution 
over  the  body  surface.  The  component  of  the  force  parallel 
to  the  direction  of  flow  is  called  the  drag  and  the  normal 
component  is  called  the  lift.  An  ideal  fluid  maintains  a 
slip  condition  at  the  solid  wall  and  does  not  give  rise  to  a 
drag  force.  However,  this  is  in  glaring  contradiction  to 
experimental  observation  where  a drag  force  is  always 
exerted  on  a body  having  a relative  motion  with  respect  to 
fluid.  Treatment  of  real  fluids  takes  into  account  a no-slip 
boundary  condition  and  the  viscosity  of  the  fluid.  Although 
the  existence  of  viscosity  had  been  considered  in  earlier 
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works,  Prandtl  (1904)  incorporated  it  in  the  momentum 
equation  by  pioneering  the  boundary  layer  analysis. 

The  discussion  on  viscosity  can  be  best  understood  by 
drawing  an  analogy  with  other  resistive  properties  of 
material.  All  real  substances  offer  certain  resistance  to 
deformation.  In  the  case  of  a solid,  if  a shear  stress  is 
applied,  then  it  is  proportional  to  the  strain  and  the 
constant  of  proportionality  (the  resistive  factor)  is  called 
the  modulus  of  torsion.  In  fluid,  the  shear  stress  is  found 
to  be  proportional  to  the  rate  of  strain  and  the 
proportionality  constant,  is  called  viscosity.  There  are 
various  concepts  to  define  and  describe  viscosity  of  fluids, 
but  the  fundamental  one  as  described  by  Newton  is  discussed. 
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Figure  1-2  Velocity  distribution  between 

parallel  plates 

To  define  viscosity  as  obtained  from  Newton's  law,  we 
can  visualize  the  motion  of  a fluid  between  two  very  long 
parallel  plates  one  of  which  is  at  rest  while  the  other 
moves  with  a constant  velocity  parallel  to  itself  as  shown 
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in  Figure  1-2.  From  experiment  it  is  observed  that  the  fluid 
adheres  to  the  walls  i.e.  u (@  y=0)  = 0 and  u (§  y=h)  = U. 

If  the  velocity  distribution  is  linear,  then  u(y)  = (y/h)U. 
In  order  to  support  the  motion,  a tangential  force  needs  to 
be  applied  to  the  upper  plate  and  this  will  be  in 
equilibrium  with  the  frictional  forces  in  the  fluid.  From 
experiment  this  force,  taken  per  unit  area,  is  proportional 
to  velocity  U and  varies  inversely  with  the  distance  h. 

Hence  the  frictional  force  per  unit  area,  denoted  by  t , is 
also  proportional  to  U/h  or  generally  du/dy.  The 
proportionality  factor  between  t and  du/dy  is  denoted  by 
and  is  a property  of  the  fluid. 


(1-1) 


The  quantity  '/j'  is  called  the  viscosity  of  the  fluid  and 
Equation  1-1  is  known  as  Newton's  Law  of  Friction. 

The  physical  concept  behind  this  definition  is  that  the 
transfer  of  momentvim  across  any  plane  in  the  fluid  is  due  to 
an  exchange  of  molecules  between  fluid  layers  on  either  side 
of  the  plane.  As  a molecule  crosses  from  one  side  with  a 
certain  velocity  to  the  other  side  of  the  plane,  it  changes 
its  velocity  if  collision  occurs.  A change  of  momentum  on 
either  side  of  the  plane  gives  rise  to  a force  parallel  to 
the  plane.  The  force  is  proportional  to  the  difference  in 
the  flow  velocity,  given  by  the  velocity  with  which  the 
molecule  started  and  the  velocity  which  it  attained  after 
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collision.  The  fluid  offers  resistance  to  this  shearing 
force  by  virtue  of  its  transport  property  called  viscosity. 
It  should  be  noted  that  the  above  description  is  a simple 
case  in  one  direction.  The  generalized  case  is  offered  by 
the  Stoke 's  law  of  friction.  The  derivation  of  the  stress 

system  and  the  rate  of  strain  for  an  infinitesimal  fluid 

* 

element  has  been  dealt  with  in  detail  by  Schlichting  (1979). 
The  relationship  as  obtained  from  Newton's  law  of  friction 
can  be  easily  deduced  from  the  general  case  as  derived  in 
the  Stokes'  law  of  friction.  One  should,  however,  note  that 
the  fluids  obeying  Stokes'  law  are  isotropic  and  Newtonian. 

A fluid  is  said  to  be  isotropic  when  the  relation  between 
the  components  of  the  stress  and  those  of  the  rate  of  strain 
is  the  same  in  all  directions  and  Newtonian  if  the  relation 
is  linear.  Fortunately,  all  gases  and  many  liquids  of 
interest  belong  to  this  class. 

It  can  be  concluded  that  viscosity  is  a measure  of  the 
fluid's  tendency  to  dissipate  energy  when  disturbed  from 
equilibrium  by  the  imposition  of  a flow  field,  which 
distorts  the  fluid.  It  appears  as  a proportionality  constant 
when  the  shear  stress  is  related  to  the  velocity  gradient  in 
a linear  way.  While  defining  viscosity  on  the  basis  of 
molecular  exchange,  the  existence  of  transport  of  any  mass 
in  a direction  normal  to  the  velocity  on  a scale  larger  than 
a molecular  one  is  excluded.  Thus,  the  definition  is 
strictly  restricted  to  laminar  flow  of  fluids. 
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1.3  Viscosity  as  an  Important  Transport  Property 
In  the  above  discussion  it  is  observed  that  in  order  to 
have  the  shear  stress,  there  must  be  a velocity  gradient  in 
the  fluid.  The  momentiam  exchange  across  this  gradient  gives 
rise  to  the  definition  of  viscosity.  This  correlation  can  be 
associated  with  other  phenomenological  relationships  such  as 
those  pertaining  to  heat  and  mass  transfer.  In  the  latter 
cases  there  are  coefficients  that  relate  the  flux  of  a 
quantity  to  its  gradient.  In  the  case  of  heat,  the 
coefficient  is  called  thermal  conductivity  while  it  is 
termed  mass  diffusivity  in  the  situation  of  mass  transfer. 
These  properties  of  fluids  (viscosity,  thermal  conductivity 
and  mass  diffusion  coefficient)  are  imbedded  in  the  rate 
equation  describing  the  transport  of  a quantity  relative  to 
the  medium  through  which  it  is  passing. 

The  existence  of  viscosity  and  its  incorporation  into 
the  solution  of  the  fluid  momentxim  equation  allowes  us  to 
treat  a real  fluid  different  from  an  ideal  one.  Although 
ideal  fluid  theory,  which  neglects  two  important  properties 
namely  compressibility  and  viscosity,  provides  excellent 
results  for  air  and  water,  it  is  unable  to  resolve  the  so 
called  d'Alembert's  paradox.  The  paradox  arises  since  ideal 
fluid  can  exert  only  normal  forces  and  this  cannot  explain 
the  drag  force  on  a body.  Prandtl  solved  this  disturbingly 
mysterious  subject  of  aerodynamics  with  his  work  on  boundary 
layer  theory  from  1904  onward.  Omission  of  viscosity  from 
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the  momentum  equation  near  the  wall  is  no  longer  permissible 
(boundary  layer  equation)  since  the  fluid  viscosity  now 
affords  the  calculation  of  drag  force,  which  is  essential  in 
designing  various  aerodynamical  bodies,  fluid  machineries, 
etc.  The  numerical  value  of  viscosity,  at  various 
thermodynamic  conditions,  has  become  increasingly  important 
in  order  to  perform  and  simulate  various  aerodynamical 
studies  and  to  calculate  momentum  in  real  fluids. 

1.4  Gas  Viscosity  Measurement  Techniques 

The  theory  and  practice  of  viscometry  is  based  upon  the 
hypothesis  postulated  by  Sir  Isaac  Newton  in  his 
"Principia" , concerning  the  magnitude  of  force  required  to 
overcome  viscous  resistance.  Following  Newton's  fundamental 
work,  the  hydrodynamics  of  the  next  century  were  concerned 
only  with  the  flow  of  ideal  or  nonviscous  fluids.  It  was  not 
until  1823  that  Navier  first  stated  the  general  equations 
for  the  motion  in  real  fluids.  Similar  equations  were 
deduced  by  Poisson  (1831)  and  by  Stokes  (1845),  who 
integrated  them  to  give  the  velocity  distribution  in  two 
cases  of  great  importance  in  connection  with  the  measurement 
of  viscosity,  namely,  flow  through  a cylindrical  tube  and 
motion  between  co-axial  cylinders. 

Meanwhile  on  the  experimental  side  both  Hagen  (1839) 
and  Poiseuille  (1846),  working  independently,  obtained  a 
relationship  between  the  flow  rate,  pressure  drop  and 
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dimension  of  capillary  tubes.  This  paved  the  path  for 
Wiedemann  (1856)  and  later  Hagenbach  (1860)  to  compute  the 
viscosity  of  fluid  flowing  through  the  tube.  In  1890  Couette 
devised  a different  method  of  measuring  viscosity  based  on  a 
system  of  two  concentric  cylinders;  viscosities  deduced  in 
this  way  gave  figures  identical  with  those  of  Poiseuille's. 
Capillary  tube  and  concentric  cylinder  viscometers  are  still 
the  primary  instruments  employed  these  days  for  viscosity 
measurement.  In  fact,  measurement  by  the  oscillating  disk 
viscometer  (a  version  derived  from  the  concentric  cylinder 
viscometer)  gives  the  standard  values  of  viscosity  for  most 
of  the  common  gases  at  room  temperature.  As  a result,  a 
short  description  of  these  viscometers  is  included  below. 

1.3.1.  Capillary  Viscometer 

The  capillary  or  transpiration  viscometer  (Figure  1-3) 
has  been  used  as  a primary  instrument  for  viscosity 
measurements  for  well  over  a century.  Although  great  care 
had  been  taken  in  the  design  of  such  an  instrument,  the 
evaluation  of  the  viscosity  has  sometimes  been  carried  out 
without  the  aid  of  adequate  working  equations,  leading  to 
inaccurate  experimental  results,  despite  the  high  precision 
of  measurements.  Experimental  measurements  are  carried  out 
by  measuring  the  pressure  drop  and  flow  across  a fully 
developed  flow  condition  in  a capillary  tube. 
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Figure  1-3  Capillary  tube  viscometer 

For  an  incompressible  fluid,  the  working  equation  for  a 
capillary  viscometer  is  given  as 

AP^8(L+na)  ^ d_2) 

where  p = pressure  drop  across  length  L 
m = mass  flow  rate  of  the  fluid 
a = radius  of  the  capillary  tube 
mQ  = kinetic  energy  correction  term  = 1.14 
n = modified  kinetic  energy  correction  term  including 
reservoir  inlet  shape  = 0.69  ± 0.004. 
p = density  of  the  fluid. 


The  above  equation  holds  good  for  0.5  £ Re  £ 100,  and  the 
entrance  to  the  capillary  is  sharp  and  square.  Popular  gas 
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viscometers  employing  capillary  tube  flow  are  Schultze's 
viscometer,  Rapp's  apparatus,  Rankine's  viscometer  and 
Edward's  constant  volvime  viscometer. 

1.3.2  Oscillatina-Bodv  Viscometer 

These  viscometers  constitute  a class  where  one  body 
oscillates  within  a hollow  body  which  is  at  rest  (Figure  1- 
4).  The  concentric  cylinder  viscometer  belongs  to  this 
family  of  viscosity  measuring  devices.  It  consists  of  an  > 
axially  symmetric  body  which  performs  torsional  oscillations 
in  a fluid.  The  system  is  suspended  from  an  elastic  wire  and 
the  changes  in  the  frequency,  caused  by  the  density  and 
viscosity  of  the  fluid,  are  measured  and  correlated  to  the 
viscosity  of  the  fluid. 


Figure  1-4  Oscillating  body  viscometer 

(fluid  external  and  internal  to  the 
body) 
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The  theory  of  all  oscillating  body  viscometers  can  be 
formulated  by  a single  mathematical  schema  (Kestin,  1988). 

JQ^[  +2A^”Jlj.+(l+Ao)g(T)  ]=M(t)  (1-3) 

Here  a(r)  is  the  angular  displacement  of  the  body  from 
equilibrium,  w = (2.7t/Tq)  is  the  natural  angular  frequency 
of  oscillation  in  a vacuum,  Tq  is  the  corresponding  period 
of  oscillation,  t = (ot  is  the  dimensionless  time,  I is  the 
moment  of  inertia  of  the  solid  body,  is  the  spring 

constant,  M(t)  is  the  frictional  moment  exerted  by  the  fluid 
and  ^ denotes  the  logarithmic  decrement  in  a vacuum  (for 
M(t)  = 0) . 

Apart  from  the  primary  instruments,  there  are  other 
devices  that  do  not  satisfy  the  accuracy  of  primary 
instrxments  but  offer  an  opportunity  to  study  a range  of 
fluid  states  not  available  to  the  other  techniques.  The 
Torsional  quartz-crystal  viscometer  and  the  falling  body 
viscometer  are  two  important  members  of  the  secondary  class. 

1 . 4 An  Introduction  to  Drag  Viscometer 
Viscosity  measurement  for  most  viscometers  is  computed 
by  ascertaining  some  parameter  which  can  be  precisely 
obtained  experimentally.  These  parameters  and  the  viscosity 
of  the  fluid  are  correlated  through  some  accurate 
expressions  which  are  derived  analytically  or  numerically. 
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Viscosity  of  a fluid  manifests  its  existence  by  providing 
drag  on  the  body  over  which  the  fluid  flows.  However,  the 
primary  instruments  do  not  measure  the  drag  force  but 
instead  some  secondary  phenomenon  associated  with  it. 
Viscosity  measurement  by  utilizing  the  drag  force  has  been 
limited  to  the  falling  sphere  viscometer  and  liquids.  The 
Reynolds  niimber  for  such  an  exercise  is  in  the  creeping 
region,  to  exclude  boundary  layer  separation  and  to  derive 
an  exact  analytical  expression.  Even  then  some  difficulty 
arises  as  the  body  may  spin  and  not  fall  straight  and 
deviate  from  the  analytical  conditions.  Moreover,  to  obtain 
this  small  Reynolds  number  for  gases,  one  has  to  employ 
micron-sized  spheres.  With  this  size  constraint  even  the 
drag  force  is  in  the  microgram  range.  Cosequently,  many 
attempts  to  measure  gas  viscosity  through  direct  drag 
measurement  were  unsuccessful.  With  the  advent  of  the  modern 
microbalance,  this  approach  becomes  feasible  and  has  been 
employed  in  this  work.  The  gas  is  allowed  to  flow  over  a 
sphere  hung  at  the  axis  of  a cylinder.  Drag  force  is 
measured  by  an  ultra-  sensitive  thermal  microbalance.  A 
relationship  has  been  obtained  to  compute  the  viscosity  from 
the  drag  force  through  dimensional  analysis. 

The  set  up  of  the  drag  viscometry  system,  has  several 
advantages  when  it  comes  to  measuring  corrosive  or  extremely 
high  temperature  gases.  The  primary  instriiments  have  been 
very  successful  up  to  a temperature  of  2000QC;  beyond  this 
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unascertained  discrepancies  take  place.  It  is  mainly  because 
the  measuring  probes  themselves  are  exposed  to  harsh 
conditions.  Apart  from  this  major  shortcoming,  in  capillary 
tube  viscometers  there  are  other  serious  problems  like  the 
exact  measurement  of  bore  diameter  due  to  thermal  expansion, 
maintenance  of  an  even  temperature  throughout  the  length  of 
the  capillary  tube,  ideal  behavior  of  gases  for  the 
governing  equations  to  hold  good  and  in  the  case  of  coiled 
tubes,  the  maintenance  of  Dean  nvimber  below  the  critical 
value  of  10  for  insignificant  secondary  flow.  In  all 
oscillating  body  viscometers  some  of  the  associated  high 
temperature  problems  include  selection  of  a suspension 
strand  which  has  long-time  elastic  stability  and  low 
internal  damping  rate,  oscillation  time  if  measured  through 
an  optical  telescope  must  take  effects  of  refraction  into 
account,  and  accounting  for  thermal  expansion  in  the 
governing  equation.  As  we  shall  see  in  the  case  of  a drag 
based  viscometer,  we  do  not  need  to  obtain  an  exact 
analytical  solution.  Results  from  dimensional  analysis  help 
in  obtaining  the  value  of  viscosity  quite  accurately.  The 
main  measuring  device,  the  thermal  microbalance,  is  never 
exposed  to  the  gas  and  is  shielded  from  the  unknown 
environment.  The  suspended  sphere  is  subjected  to  the  high 
temperature  and  a temperature  gradient  exists  in  the  thermal 
boundary  layer  close  it.  As  a result,  the  gas  has  a much 
lower  temperature  after  it  has  flown  over  the  sphere. 


CHAPTER  2 

SURVEY  OF  PREVIOUS  EFFORTS  FOR  DRAG  MEASUREMENT  ON  SPHERE 

2 . 1 Drag  on  an  Unbounded  Sphere 
The  drag  force  on  a body,  exerted  by  the  fluid  flowing 
over  it,  is  obtained  by  the  solution  of  the  Navier-Stokes ' 
equation.  However  an  exact  result  is  not  obtained  in  most  of 
the  cases.  Analytical  treatment  of  the  Navier-Stokes 
equation  for  fluid  flow  is  amenable  at  very  low  flow  rates. 
At  higher  flow  rates  experimental  results  have  been  obtained 
for  the  drag  force  and  analysis  is  done  through  numerical 
procedures.  Since  a host  of  physical  phenomena  occur  for  the 
flow  over  a sphere  at  different  Reynolds  niimbers,  empirical 
correlations  differ  for  different  flow  regimes.  Description 
of  the  flow  field  over  a sphere  has  been  provided  in  several 
text  books.  Broadly  speaking  the  flow  over  a sphere  remains 
unseparated  until  a Reynolds  number  of  around  20.  Beyond 
this  there  is  onset  of  separation  and  the  resulting  wake 
remains  steady  until  a Reynolds  number  of  130.  After  this 
discrete  pockets  of  vortices  begin  to  be  shed  and  the  wake 
looks  unstable.  For  Reynolds  numbers  from  about  20  to  400, 
the  separation  angle  in  degrees  from  the  front  stagnation 
point  is  given  by  0g  = 180  - 42.5[ln(Re/20]°-^®^.  At  Reynolds 
numbers  beyond  400  the  unsteady  wake  also  becomes 
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asymmetric.  Throughout  the  separation  angle  moves  from 
behind  the  sphere  in  a forward  direction;  i.e.,  it  moves 
against  the  flow  direction.  But  at  a Reynolds  number  beyond 
2x10^  a drastic  change  in  flow  pattern  takes  place  and  the 
separation  begins  to  move  aft. 

There  is  no  single  relation  that  is  able  to  describe 
these  different  flow  behaviours.  For  the  case  of  a sphere 
G.G. Stoke  (1845)  was  the  first  to  give  a solution  for  steady 
creeping  flow  past  a rigid  sphere.  The  analytical  treatment 
was  possible  by  neglecting  the  convective  terms  of  the 
momentiim  equation.  C.W.Oseen  (1910)  realized  that  the  Stokes 
equation  for  creeping  flow  and  finite  Reynolds  number  is 
valid  for  distances  less  than  a/Re,  where  a is  the  radius  of 
the  sphere  and  Re  the  Reynolds  number.  He  incorporated  the 
neglected  inertial  terms  in  a linearized  form.  Several 
series  solutions  were  obtained  for  this  Oseen  equation  by 
including  from  1 to  24  terms.  Following  another  line  of 
thought  I.  Proudman  and  J.  R.  A.  Pearson  (1957) used 
perturbation  techniques  to  solve  the  creeping  flow  over  a 
sphere.  But  all  these  analytical  exercises  were  limited  to  a 
Reynolds  number  of  less  than  unity. 

The  boundary  layer  theory,  which  is  an  approximation  of 
the  Navier-Stokes  equation  obtained  by  neglecting  the 
streamwise  diffusion,  has  gained  partial  success  in 
predicting  fluid  velocity.  But  it  is  applicable  at  higher 
Reynolds  nximbers  (Re>3000)  when  the  boundary  layer  thickness 
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is  thin.  There  are  severe  limitations  to  this  approach 
although  it  does  provide  a meaningful  way  to  explain  the 
fluid  flow  phenomenon.  The  major  shortcomings  are  that:  1) 
the  outer  boundary,  which  is  treated  as  potential  flow,  is 
true  to  the  assumption  until  30°  from  the  front  stagnation 
point;  2)  once  the  flow  separates,  there  is  a convincing 
feedback  by  the  downstream  on  to  the  oncoming  flow;  and  3) 
the  boundary  layer  may  not  be  negligible  compared  to  body 
size  at  low  flow  rates.  The  latter  effect  has  been  accounted 
for  by  Gluckman  (1971)  through  incorporation  of  a 
displacement  thickness.  Inspite  of  all  this,  not  much 
headway  has  been  made  in  this  field  for  separated  flows. 

Numerical  solution  to  the  Navier-Stokes  equation  has 
been  carried  out  to  get  drag  measurement  at  higher  flow 
rates  when  separation  occurs.  Due  to  computing  difficulties 
encountered  in  the  earlier  days,  the  treatment  was  limited 
to  a Reynolds  nvimber  for  which  asymmetry  had  not  set  in. 
Several  approaches  to  solve  the  Navier-Stokes'  equation  by 
the  stream  function-vorticity  method  were  undertaken.  Jenson 
(1959)  and  LeClair  et  al.  (1970)  used  a steady  state 
equation  while  Rimon  and  Cheng  (1969)  and  Rafique  (1971) 
used  a time  dependant  momentum  equation.  Dennis  and  Walker 
(1971)  expanded  stream  function  and  vorticity  as  a series  of 
Legendres  functions  and  their  substitution  into  the 
governing  equation  resulted  in  a series  of  ordinary 
differential  equations.  However,  this  approach  is  not 
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attractive  at  higher  Reynolds  numbers  since  the  number  of 
terms  to  be  included  becomes  prohibitively  large. 

Thus  far  the  experimental  study  has  unequivocally 
provided  reliable  values.  The  drag  forces  experimentally 
obtained  have  been  presented  in  the  form  of  a plot  of 
coefficient  of  drag  (Cp)  vs.  Reynolds  number  (Re)  in  several 
fluid  mechanics  books.  In  this  study,  the  fluid  situation  is 
limited  to  laminar  flows.  Some  of  the  important  drag 
correlations  in  this  regime  are  (Clift  et  al.,  1978) 

Cr>=^—+—  ; Re<0.01  (2-1) 

° 16  Re  ' 

Cn  = .^[  1+0. 1315Re(°*®2‘°'°^*'h  ; 0.01^Re£20  (2-2) 

^ Re'- 

[1+0. 1935Re°*®^°®]  ; 20^Re^260  (2-3) 

Rq 

logioC^=l.  6435-1. 124w+0.1558w2  ; 260^J?e^l500  (2-4) 

where  w=log2®*  Other  correlations  for  the  drag  forces  have 
been  judiciously  compiled  by  R. Clift  et  al.  (1978),  and  it 
was  deemed  unnecessary  to  duplicate  the  information  here. 

In  conclusion,  for  free  flow  over  a sphere,  only 
experimental  and  numerical  techniques  can  be  used  to  provide 
accurate  results  at  a high  Reynolds  number.  In  the 
experimental  field  the  state-of-the-art  instruments  and 
measuring  techniques  will  continue  to  provide  accurate 
values.  On  the  numerical  side,  both  the  finite  difference 
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method  and  finite  element  method  are  being  used  as 
proficient  techniques  to  gain  success  in  this  area. 

2.2  Drag  on  a Bounded  Sphere 
The  flow  of  fluid  over  a sphere  confined  by  an  external 
wall,  is  further  complicated  by  the  wall  effects.  As 
compared  to  an  unbounded  sphere,  little  work  has  been 
reported  in  this  field.  Theoretical  work  has  been  done  for 
the  creeping  flow  range.  Coutanceau  (1972)  has  provided  some 
flow  visualization  results  and  Johansson  (1974)  worked  out  a 
numerical  approach.  These  studies  have  been  very  limited  in 
their  function  because  of  the  various  difficulties 
encountered. 

In  general,  it  has  been  found  that  the  flow  situation 
is  governed  by  the  size  of  the  sphere,  the  cylinder  size, 
the  position  and  surface  roughness  of  the  sphere  as  well  as 
the  velocity  profile  of  the  oncoming  flow  (Figure  2-1).  The 
effect  of  sphere  and  cylinder  diameter  size  has  been 
directly  related  in  various  correlations  and  this  is  noted 
later.  The  position  of  the  sphere  within  the  tube  gives  rise 
to  an  additional  pressure  drop  within  the  latter.  For 
standard,  fully  developed  flow  a smooth  sphere,  the  drag 
forces  on  the  sphere  and  the  wall  as  a function  of  the 
sphere  position  has  been  looked  into.  Without  going  into 
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Figure  2-1  Sphere  position  and  fluid  motion  in  cylinder 

mathematical  details,  it  has  been  shown  (Happel  and  Brenner, 
1973)  that  to  obtain  a first  approximation  to  the  increased 
pressure  drop  arising  from  the  presence  of  a particle  in  a 
bounded  flow,  the  momentum  theorem  gives 

F>(^-1)F  (2-5) 

where  is  the  extra  shearing  force  on  the  cylinder  wall 
due  to  the  presence  of  the  sphere,  Vjj,  is  the  mean  velocity 
with  respect  to  the  wall,  Vg  is  the  approach  velocity  to  the 
sphere  and  F is  the  drag  force  on  the  sphere.  For  a circular 
cylinder  of  radius  R 
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^0 

Vm 


(2-6) 


Thus  from  Equation  2-5  we  find  that  if  one  moves  from  the 
cylinder  axis  (b/R  =0)  towards  the  wall  (b/R  =1),  the 
direction  of  the  extra  force  changes  from  being  parallel  to 
the  drag  to  being  oppositely  directed.  As  a special  case, 
when  b/R  = 0.707,  the  force  no  longer  exists.  For  the  case 
of  determining  the  effects  of  surface  roughness  of  the 
sphere,  not  much  has  been  done  at  low  to  moderate  flow 
rates.  Achenbach  (1974)  has  determined  over  the  Reynolds 
nvimber  range  (based  on  the  minimum  cross  sectional  area  at 
the  equator  of  the  sphere)  of  3xlO^<Re<2xlO®  the  effects  of 
surface  roughness.  He  showed  that  increasing  the  roughness 
parameter  increased  the  critical  drag  coefficient  to  a 
maximum  value  of  0.4.  However  the  Strouhal  nuinber  for  the 
vortex  shedding  remained  unaffected  and  was  equal  its  value 
for  a smooth  sphere. 

Most  of  the  work  has  been  performed  by  experimentation. 
The  results  obtained  can  be  divided  on  the  basis  that  the 
sphere  is  falling  in  the  fluid  or  it  is  held  fixed  against 
an  oncoming  fluid.  In  either  case,  studies  have  been 
performed  to  take  into  account  the  wall  effect.  A wall 
correction  factor  'K'  is  normally  used  which  gives  the  ratio 
of  the  drag  force  on  the  sphere  in  a bounded  condition  to 
that  when  the  sphere  is  in  infinite  fluid  flow. 
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2 ♦ 2 « 1 Sphere  Sedimenting  in  a Cylindrical  Tube 

Much  work  has  been  done  for  the  case  of  a sphere 
sedimenting  in  a cylindrical  tube.  It  was  started  as  early 
as  1687  by  Sir  Isaac  Newton  and  has  received  attention  since 
then.  As  in  the  case  of  an  unbounded  sphere,  there  is 
difficulty  in  deriving  a single  expression  for  the  fluid 
flow  phenomenon.  The  interference  effect  of  the  cylindrical 
walls  on  the  falling  sphere  is  a function  of  Reynolds  nximber 
and  it  decreases  with  increasing  Reynolds  number. 


Table  2-1  Correlations  for  sphere  sedimenting  in  a tube 

at  low  flow  rates  (Re<l) 


1 Author 

RELATIONSHIP  FOR  K (Cu/Ci3„« 

1 Ladenburg  (1907) 

(1+2. IB);  B<0.1 

1 Faxen  ( 1921) 

1/ ( 1-2 . 104B+2 . 09B^-0 . 95B^ ) ; B<0 . 2 

1 Francis  ( 1933 ) 

(1-0.475B)'^/(1-B)'^;  B<0.9 

1 Haberman  (1958) 

( 1-0 . 75857B® ) / ( 1-2 . 105B+2 . 0865B^- 
1.7068B®+0.72603B®) ; B<0.9 

Some  of  the  important  relationships  for  the  wall 
correction  factor  ”K"  as  a function  of  "B”  (diameter  ratio 
of  the  sphere  and  the  cylinder)  for  a sedimenting  sphere  on 
the  axis  of  the  cylinder, are  shown  in  Table  2-1.  The 
variation  in  the  coefficient  of  drag  with  Reynolds  number 
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has  been  plotted  in  Figure  2-2.  For  higher  flow  rates 
(Re>1000)  some  of  the  published  values  of  "K"  are  provided 
in  Table  2-2.  The  relationships  in  Table  2-1  were  obtained 
by  the  analytical  solution  of  the  fluid  flow  equation  at 
creeping  range  while  results  in  Table  2-2  were  through 
experimental  procedure. 

Table  2-2  Correlations  for  sphere  sedimenting  in  a tube 

at  high  flow  rates 


AUTHOR 

RELATIONSHIP  FOR  K (Cd/Cd<„) 

Newton  (1687) 

(l-fi2) (1-b2/2)0-5 

Munroe  (1888) 

I-J3I.5 

Lunnon  (1928) 

1-0. 23B  ; B<0.3 

Mott  (1951) 

1/(1+16B'^) ; 0.5<B<0.7 

1/(1+AB2);  0.2<B<0.5  and  1.8<A<3.2 

2.2.2  Sphere  Fixed  in  a Cylindrical  Tube 

For  the  case  of  a sphere  held  fixed  at  the  axis  of  the 
cylinder  and  against  an  oncoming  flow,  different 
correlations  are  available  from  the  literature.  In  order  to 
understand  the  phenomenon  of  fluidization,  sedimentation  and 
flow  through  fixed  beds,  researchers  came  up  with 
relationships  for  the  drag  force  on  an  unbounded  and  bounded 
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sphere.  The  important  ones  for  creeping  flows  are  give  in 
Table  2-3. 


Table  2-3  Correlations  for  sphere  fixed  in  a tube 

at  low  flow  rates  (Re<l) 


AUTHOR 

RELATIONSHIP  FOR  K 1 

Wakiya  (1957) 

(1-0.1667B2)/(1-2.105B  +2.09B®  1 

l.llB®)  1 

Brenner  and  Happel 

(1-0.1667B2)/( 1-2. 105B) 

(1958) 

Haberman  and  Sayre 

(1-0. 1667B2-0.20217B®)/( 1-2. 105B+ 

(1958) 

2.0865B®-1.7068B®+0.72603B®) 

For  higher  Reynolds  numbers  the  notable  relationships 
for  ••K”  are 

a)  Fayon  (1960)  K=[  ( 1-0 . 1667fi2 )/( 1-2  . lOSB+2 . 08713^ ) 

+ {(Cj^/Cg)-!}} 

where  and  Cg  are  the  actual  coefficient  of  drag  and  the 
Stokes  drag  coefficient,  respectively.  The  relationship  is 
valid  for  Reynolds  numbers  less  than  40  and  J3<0.3. 

b)  Clift  et  al.  (1978)  K=l/(  1-1. 6B^*®)  ; for  B<0.6  and 

100<Re<1000 

c)  Achenbach  (1972)  K=(  1+1 . 4B‘*- ^ ) / ( l-B^  )2  ; for  Re  > 10®. 
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A sphere  may  also  exist  in  a cylinder  without  being 
sedimented  or  fixed  by  a support;  it  may  be  freely  suspended. 
However,  these  freely  suspended  spheres  generally  do  not 
remain  fixed  in  space  with  respect  to  the  containing  tube. 
They  may  spin  about  their  own  axis  and  simultaneously  rotate 
about  the  axis  of  the  tube.  Such  dynamic  effects  are  extremely 
difficult  to  handle  analytically  and  cannot  be  predicted  by 
dimensional  analysis  as  a function  of  geometry  although  the 
phenomenon  is  itself  related  to  the  geometry.  For  these  freely 
suspended  spheres  it  has  been  found  that  for  J3=0 . 2 , the 
spheres  tend  to  stay  to  one  side  of  the  center  of  the  tube  and 
cling  to  the  wall.  Larger  diameter  spheres  spin  about  their 
own  axis  and  about  the  axis  of  the  tube.  Also,  in  spite  of  a 
fundamental  difference  in  flow  pattern  past  a fixed  and  freely 
suspended  sphere,  the  drag  coefficient  is  approximately  the 
same  for  all  practical  purposes  (Round  et  al.  1972).  However, 
experiments  conducted  this  way  to  correlate  drag  coefficient 
with  Reynolds  number  for  fixed  spheres  have  not  been 
considered  in  this  study. 
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Reynolds  Nvunber  (Re) 


Figure  , 2-2  Plot  of  Cq  vs. 

tube 


Re  for  sedimenting  sphere  in  a 


Coefficient  of  Drag  (Cd) 
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Figure  2 3 Plot  of  vs  Re  for  fixed  sphere  in  a tube 


CHAPTER  3 

EXPERIMENTAL  PROCEDURE 


3 . 1 Considerations  for  Experimentation 
For  the  results  of  any  good  experimentation  to  convey  a 
meaningful  information,  some  basic  guidelines  have  to  be 
adhered  to.  These  guidelines  vary  with  the  nature  of 
experiments.  The  current  study  has  the  objective  of 
measuring  the  drag  force  over  a bounded  and  suspended 
sphere . The  measurements  have  then  to  be  extended  for  the 
determination  of  viscosity.  In  general  for  primary 
viscometry,  the  experimental  methods  selected  should  satisfy 
the  following  conditions  (Kestin  and  Wakeham,  1988): 

a)  Maintenance  of  a near-equilibrium  state  by  use  of  small 
shear  rates. 

b)  The  preservation  of  hydrodynamic  stability  or,  at  least, 
the  reduction  of  the  effects  of  unavoidable  instabilities  to 
a negligible  level. 

c)  The  availability  of  a faithful  mathematical  model  of  the 
process  and  of  an  accurate  solution  in  the  form  of  an 
adequate  working  equation  and  a complete  set  of  corrections. 

d)  High  resolution  in  the  detection  of  the  major  effect  and 
operation  at  its  lowest  value  consistent  with  the  desired 
precision  and  the  sensitivity  of  the  working  equation. 
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Compliance  to  these  conditions  usually  helps  in 
obtaining  reliable  values  of  viscosity  for  fluids. 

3 ♦ 2 Description  and  Working  Principle  of  Apparatus 

The  experiment  was  performed  by  setting  up  the  system 
as  shown  in  Figure  3-1.  Drag  force  measurement  was  done  by 
using  a Cahn  Instrument  manufactured  C2000  model 
microbalance.  The  thermal  microbalance  is  a very  sensitive 
weight  and  force  measurement  instrviment.  It  is  designed  for 
weights  and  forces  up  to  3.5  gram-force  and  is  sensitive  to 
changes  as  small  as  0.1  microgram.  The  balance  is  divided 
into  two  sections;  one  section  is  the  control  unit  where  all 
controls  and  outputs  are  contained  and  the  other  section  is 
the  weighing  unit  which  detects  the  actual  weight  or  force . 
The  weighing  unit  may  be  operated  in  vacuum,  flowing  gas  or 
atmospheric  conditions.  The  balance  consists  of  (1)  a 
balance  beam  mounted  to,  supported  by  and  pivoting  about  the 
center  of  a taut  ribbon;  (2)  a torque  motor  coil  located  in 
a permanent  magnetic  field  and  also  mounted  to  the  taut 
ribbon;  (3)  sample  suspension  fixtures;  (4)  a beam  position 
sensor  system;  and  (5)  controls,  circuitry  and  indicators 
(Figure  3-2).  Weights  or  forces  to  be  measured  are  applied 
to  the  sample  (left)  side  of  the  beam  which  produces  a force 
about  the  axis  of  rotation.  An  electric  current  flowing  in 
the  torque  motor  also  produces  a force  about  the  same  axis 
which  is  equal  and  opposite  to  the  force  from  the  beam  if 
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Figure  3-1  Experimental  set-up  for  drag  measurement 
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Figure  3-2  Microbalance  beam  assembly 
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the  beam  is  at  the  reference  beam  position.  This  reference 
position  is  detected  by  the  beam  position  sensing  system.  A 
greater  force  on  the  beam  will  require  a greater  opposite 
force  from  the  torque  motor  to  keep  the  beam  at  the 
reference  position.  Therefore,  the  current  necessary  to 
produce  the  required  torque  motor  force  is  a direct  measure 
of  the  force  on  the  beam.  The  process  of  calibration  allows 
this  current  to  be  measured  in  units  of  mass  (grams)  or 
force.  The  sensitivity  and  readability  is  0.1  micrograms; 
repeatability  is  10“^  of  the  total  load  (max.  0.2 
microgram);  accuracy  is  ±0.1%  of  the  recorder  range  plus  1.5 
X 10”^  of  the  suppression  range  for  absolute  weighing. 

The  gas  flow  meter  used  was  an  Omega  Engineering  Inc. 
FMA  800  series  mass  flow  meter  with  a range  of  0-50  Lt/min. 
Like  most  mass  flow  sensors,  it  measures  mass  flow  directly 
without  the  need  for  temperature  or  pressure  corrections, 
and  the  controllers  can  automatically  control  flowrate 
regardless  of  variation  in  line  pressure  (over  a limited 
pressure  range ) . The  measuring  unit  features  a corrosion 
resistant  stainless  steel  and  Viton  construction  and  a 50 
micron  inlet  screen  to  prevent  contamination  of  the  sensor. 
The  flow  body  has  a leak  less  than  10  atm-cc/sec  of  heJJijaBU — 
The  FMA  800  series  mass  flow  sensor  operates  on  the 
principle  of  heat  transfer,  by  sensing  the  change  in 
temperature  (AT)  along  a heated  section  of  a capillary  tube. 
Part  of  the  total  flow  is  forced  through  the  capillary  by 
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means  of  a pressure  differential  of  approximately  0.45  psi. 
The  required  maximum  flow  rate  determines  the  number  of 
discs;  these  are  bolted  together  and  constitute  the  laminar 
flow  device.  The  design  of  the  discs  is  such  that  flow 
conditions  in  both  the  capillary  and  laminar  flow  device  are 
comparable,  thereby  resulting  in  proportional  flow  rates 
through  them.  The  At  sensed  by  the  upstream  and  downstream 
temperature  sensors  on  the  capillary  depends  on  the  amount 
of  heat  absorbed  by  the  mass  of  the  gas  flow.  The  mass  of 
the  gas  being  virtually  independent  of  the  temperature  and 
pressure  makes  the  instrument  independent  of  temperature 
and  pressure  variations  i.e.  a mass  flow  meter. 

3.3  Experimentation 

During  experimental  measurements,  consideration  was 
given  to  the  tube  size,  flow  rate  and  sphere  diameter  and 
its  weight.  Teflon  balls  of  different  diameter  (13=  0.35, 
0.42,  0.49,  0.56)  were  centrally  hung  into  a long 
cylindrical  tube  of  diameter  22.45  mm.  The  spheres  were 
suspended  by  means  of  a fine  nylon  string  (d=0.07  mm)  from 
loop  B of  the  microbalance  (Figure  3-2).  No  special  attempt 
to  determine  the  sphericity  of  the  teflon  balls  was  made. 
Drag  measurement  using  the  C2000  microbalance  was  dictated 
by  maximum  mass  restriction  of  3.5  grams  and  this  was  a 
guide  line  for  choosing  the  light  weight  and  smooth  teflon 
balls.  Operating  limits  were  calculated  approximately  to 
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insure  that  the  tube  Reynolds  number  was  small  enough  to 
avoid  either  incipient  transition  to  turbulence  or  unduly 
long  entrance  length.  A rough  estimate  for  the  entrance 
length  Lg  in  tubes  can  be  obtained  from  the  relationship 
Lg/D  = 0.06  Re;  where  D is  the  tube  diameter  and  Re  is  the 
Reynolds  number  based  on  the  mean  flow  velocity.  In  all  the 
measurements  reported  here,  this  criterion  was  satisfied. 

Drag  on  the  sphere  was  obtained  by  flowing  nitrogen  gas 
with  a Reynolds  number  (based  on  mean  velocity  and  tube 
diameter)  of  40  to  1000.  The  difference  in  reading,  as  given 
by  the  microbalance,  with  the  gas  flowing  and  with  no  gas 
flow  gives  us  the  total  drag  force  exerted  on  the  ball  and 
the  string.  The  drag  on  the  string  was  calculated  by  using 
different  string  lengths  for  hanging  the  same  ball  under  the 
same  flow  conditions.  The  difference  of  the  two  readings,  at 
two  different  string  lengths,  was  due  to  the  drag  on  the 
string  alone.  This  value  was  subtracted  from  the  original 
reading  to  obtain  the  drag  on  the  sphere  alone.  For  example; 
let  = Drag  on  sphere  (Dgpj^)  + Drag  on  first  string  (Dg^i) 
D2  = Drag  on  sphere  (Dgp2)  + Drag  on  second  string  (Dgt2) 
If  the  length  of  the  second  string  (L2)  is  greater  than  the 
first  (Lj^)  then  the  drag  on  the  string  per  unit  length  is 

/ ^spl  ~ ^1  ~ (^unit 

^sp2  “ ^2  “ (^unit  )^2 
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The  two  values  of  D obtained  in  this  way,  maintaining  the 
same  flow  conditions  for  the  sphere  and  tube,  should  be 
equal  if  the  approximation  is  good.  Table  3-1  provides  some 
of  the  values  obtained  to  justify  the  point. 

Table  3-1  Table  of  drag  values  obtained  by  two 

different  string  lengths 


Re 

B (d/D) 

^unit 

(N/m) 

^spl 

(N) 

^sp2 

(N) 

187.46 

0.35 

1.28419E-8 

2.71619E-8 

2.71619E-8 

429.32 

0.35 

1.39126E-7 

6.71933E-8 

6.71933E-8 

680.08 

0.35 

2. 16715E-7 

1.45222E-7 

1.45222E-7 

187.40 

0.42 

5.45800E-8 

2.96641E-8 

2.96641E-8 

434.38 

0.42 

1. 11728E-7 

1.01855E-7 

1.01855E-7 

678.37 

0.42 

1.63716E-7 

2. 14169E-7 

2. 14169E-7 

185.69 

0.49 

1.01428E-7 

3.76202E-8 

3.76196E-8 

430.96 

0.49 

1.50775E-7 

1.31952E-7 

1.31951E-7 

677.95 

0.49 

2.78858E-7 

2.73101E-7 

2.73097E-7 

194.84 

0.56 

1.11124E-7 

5.58023E-8 

5.58022E-8 

431.15 

0.56 

2. 15246E-7 

1.91915E-7 

1.91915E-7 

675.27 

0.56 

3.54523E-7 

4.05596E-7 

4.05596E-7 
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In  Table  3-1  the  Reynolds  number  is  based  on  tube  diameter 
and  mean  velocity  of  flow. 

Therefore  it  is  safe  to  assume  that  the  method  to 
obtain  value  of  drag  force  on  the  sphere  by  using  two 
different  string  lengths  is  quite  accurate.  Figure  3-3  shows 
the  coefficient  of  drag  obtained  during  this  experiment  for 
sphere  to  cylinder  diameter  ratio  (6)  of  0.5,  and  as  given 
by  other  authors. 


Figure  3-3  Coefficient  of  drag  (Co)  vs.  Reynolds  number 

(6=0.5) 
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3 . 4 Experimental  Observations 
In  order  to  conduct  a meaningful  experiment  it  is 
necessary  to  obtain  a steady  operational  point.  Some  of  the 
salient  features  which  were  observed  during  the  experimental 
runs  and  which  in  turn  helped  us  to  set  up  a good 
experimental  platform  are  briefly  discussed  here.  It  was 
noticed  that  if  the  sphere  was  hung  too  low  in  the  cylinder 
two  problems  arose:  a)  the  undeveloped  flow  caused 
oscillation  in  the  drag  reading  output  and  b)  the  drag  on 
the  string,  which  should  be  low  for  approximate  elimination, 
becomes  comparable  to  the  sphere  drag.  This  problem  is 
overcome  by  utilizing  a fully  developed  flow  condition.  In 
this  situation  the  readings  obtained  are  steady.  It  was  also 
observed  that  the  repeatability  of  the  experiment  was 
dependant  on  the  ball  diameter.  When  B (=d/D)  is  equal  to 
0.5,  experimental  results  repeatedly  matched  the  same 
values.  Although  not  much  work  was  done  to  corroborate  it, 
an  eccentric  position  of  the  sphere  seemed  to  affect  the 
readings  by  making  it  unstable.  Thus  care  should  be  taken  to 
hang  the  sphere  in  the  central  position.  Lastly,  elimination 
of  string  drag  by  using  two  different  string  lengths,  under 
the  same  set-up,  provided  excellent  results  for  calculating 
the  drag  on  the  sphere  alone.  Drag  measurements  by  this 
tactic  gave  matching  results  for  sphere  drag  with  two 
different  string  lengths  for  4 decimal  places  (in  mg  scale). 
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This  helps  us  to  accurately  plot  the  dependence  of  the 
coefficient  of  drag  on  flow  Reynold's  number. 

For  viscosity  measurements  the  following  points  need  to 
be  addressed.  Firstly,  a low  flow  rate  is  needed.  This 
assures  that  the  fluid  has  deviated  a negligible  amount  from 
its  equlibrium  condition.  Since  viscosity  is  a 
thermophysical  property,  it  is  dependant  on  the 
thermodynamic  condition.  In  order  to  measure  it,  a velocity 
gradient  needs  to  be  created.  This  dissipation  creates  a 
temperature  change,  however  small  it  may  be.  To  keep  the 
change  to  a minimum,  creeping  flows  are  normally  employed. 

Thus,  good  measurements  are  performed  in  the  presence  of 

very  small  gradients  and  the  average  of  pressure, temperature 

and  density  values  suffice  to  indicate  the  local 

thermodynamic  state.  Our  experiments  have  been  performed  at 

higher  Reynolds  number  to  basically  understand  the  physics 

at  several  flow  rates.  However,  the  accuracy  of  the 

measuring  instriiment  is  capable  of  making  very  sensitive  and 

precise  measurements  at  low  flow  rates.  Secondly,  the 

preservation  of  hydrodynamic  stability  is  attained  through 

fully  developed  flow,  using  a sphere  diameter  which  is  half 

that  of  the  cylinder  (B  = 0.5),  and  the  central  positioning 

of  the  sphere  in  the  cylinder.  Thirdly,  an  exact 

mathematical  expression  for  the  drag  force  exerted  on  a 

sphere  centrally  suspended  in  a tube  and  experiencing  a 

Poiseuille  flow  is  given  by  the  Haberman  and  Sayre  (1958)  equation 
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(l-lp2-0.20217p^) 

DRAG  - ■>— 

(1-2. 105 p+2. 0865 p^ -1.7068 p® +0.72603 p^) 

(3-1) 

The  negative  sign  indicates  that  the  drag  is  opposite  to  the 
direction  of  the  flow.  And  finally,  the  high  resolution  in 
the  detection  of  the  major  effect,  which  happens  to  be  the 
drag  force,  is  deftly  acquired  through  modern  ultra  precise 
instrximents . All  these  properties  contribute  to  making  the 
measurement  of  viscosity  through  drag  force  measurement  a 
successful  approach. 

3 . 5 Viscosity  from  Drag  Measurement 
The  total  drag  on  an  object  in  a stream  of  fluid 
consists  of  skin  friction  and  the  form  or  pressure  drag.  The 
former  is  the  integral  of  all  the  shearing  stresses  taken 
over  the  surface  of  the  body  while  the  latter  is  the 
integral  of  normal  force.  The  sum  of  the  two  is  called  the 
total  or  profile  drag.  For  free  flow  the  skin  friction  can 
be  calculated  with  some  accuracy  by  employing  boundary  layer 
theory  and  it  is  here  that  the  role  of  viscosity  also  comes 
into  the  picture.  The  form  drag  which  does  not  exist  in 
frictionless  flow,  is  due  to  the  fact  that  the  presence  of  a 
boundary  layer  modifies  the  pressure  distribution  on  the 
body  as  compared  to  an  ideal  fluid.  However,  the  computation 
of  form  drag  is  very  difficult.  Consequently,  reliable  data 
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on  total  drag  is  in  general  obtained  by  measurement.  Once 
the  total  drag  (profile  drag)  is  obtained,  the  coefficient 
of  drag  (Cq)  for  the  body  can  be  calculated  from: 

D = lpU^ACo  (3-2) 

where  D is  the  profile  drag,  p is  the  density  of  the  fluid, 

U is  the  velocity  of  the  fluid  and  A the  cross  sectional 
area  of  the  body.  From  dimensional  analysis,  we  find  that 
for  a particular  configuration  and  flow  condition  there  is  a 
unique  curve  for  Cq  and  the  flow  Reynolds  number  (Re).  The 
main  aim  is  to  accurately  plot  this  curve.  After  this  curve 
has  been  obtained,  viscosity  can  be  calculated  as  follows: 

a)  The  sample  gas  is  flown  over  the  bounded  sphere  and  the 
mean  flow  velocity  is  accurately  determined  by  a flowmeter. 
Care  should  be  taken  so  that  the  sphere  is  centrally 
positioned  in  the  tube  and  the  velocity  is  fully  developed. 

b)  From  equation  3-2  the  coefficient  of  drag  (Cq)  is 
calculated.  It  is  assximed  that  the  density  of  the  gas  is 
known  and  so  is  A,  the  cross  sectional  area  of  the  sphere. 

c)  From  the  plot  of  vs.  Re,  the  Reynolds  nximber  is 
obtained  which  gives  the  value  of  the  viscosity. 

Values  of  viscosity  were  obtained  for  oxygen  and  helium 
gas  at  room  temperature  by  this  method.  Table  3-2  shows  the 
some  of  the  errors  that  were  obtained  during  this  study. 
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Table  3-2  Errors  obtained  at  room  temperature  during 

viscosity  measurements 


Reynolds  Number 

©2  Error  ( % ) 

He  Error  (%) 

150.0 

0.77 

0.92 

200.0 

0.71 

0.88 

250.0 

0.70 

0.90 

300.0 

0.68 

0.85 

400.0 

0.65 

0.83 

500.0 

0.82 

0.93 

CHAPTER  4 
PROBLEM  ANALYSIS 


4.2  Drag  Force  on  a Cylindrically  Bounded  Sphere 

A short  description  of  this  problem  has  already  been 
provided  in  Section  2.2.2.  Since  this  happens  to  be  the 
backbone  for  comparison  of  our  results  with  different  work, 
some  more  information  is  provided  here. 

Due  to  the  complexity  of  the  general  Navier-Stoke's 
equation,  analytical  derivation  of  the  drag  force  has  been 
limited  to  the  creeping  flow  regime.  In  this  region  the 
Navier-Stokes  equation  reduces  to  a linear,  fourth  order, 
partial  differential  equation.  Initial  treatment  concerned 
uniform  flow  (Lamb  1932)  and  later  boundary  conditions  were 
imposed  by  considering  an  infinitely  long  cylindrical  tube 
(Lee  1947;  Haberman  & Sayre  1958).  Finally  solution  for 
creep  flow  about  a sphere  bounded  by  cylindrical  walls  and 
suspended  in  a Poiseuille  flow  was  undertaken.  Johanssen 
(1974)  reported  numerical  calculations  of  flow  around  a 
sphere  fixed  on  the  axis  of  a Poiseuille  flow.  Only 
solutions  for  J3  (d/D)  = 0.1  were  reported.  The  term  J3  is  the 
diameter  ratio  of  the  sphere  (d)  and  the  cylinder  (D).  Some 
experiments  were  carried  out  at  slightly  higher  Reynolds 
number.  Fayon  and  Happel  (1960)  performed  experiments  upto 
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Reynolds  numbers  of  40  while  Eichhorn  and  Small  (1964) 
extended  the  study  upto  Reynolds  numbers  of  250.  However, 
these  works  were  handicapped,  as  concluded  by  their  authors 
themselves,  in  not  measuring  the  drag  accurately  and  in 
their  inability  to  change  some  important  variables  which 
were  related  to  the  operating  characteristics  of  the 
apparatus  and  intrinsic  to  the  experimental  set-up. 

At  the  other  end  of  the  flow  spectrxim,  experiments  have 
been  performed  for  flow  past  spheres  at  very  high  Reynolds 
number  (Re  > 10^)  (McNown  & Newlin  1951;  Achenbach  1974). 

The  basic  attempt  was  to  relate  the  coefficient  of  drag  to 
the  flow  Reynolds  number  based  on  the  mean  approach  velocity 
and  B (=d/D).  The  cylindrical  wall  increases  the 
supercritical  drag  coefficient  well  above  the  value  of  0.3, 
as  arbitrarily  chosen  to  define  critical  number  for  a sphere 
in  an  unbounded  flow.  The  critical  Reynolds  number  based  on 
the  mean  approach  velocity  decreases  from  3.65  X 10^  in  an 
unbounded  fluid  to  1.05  X 10®  for  d/D=0.916. 

Taking  these  results  into  account,  Clift  et  al  (1978) 
have  tried  to  bridge  the  region  between  the  high  and  low 
flow  regimes  by  interpolation.  However,  this  curve  fitting 
method  does  not  model  any  actual  experimental  data  and  the 
results  are  open  to  doubts.  The  existing  experimental 
findings  in  the  laminar  region  (Klyachko  1934;  Young  1960) 
offer  results  which  are  not  in  agreement.  This  is  primarily 
due  to  the  varied  experimental  approach,  the  manner  in  which 
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inlet  flow  conditions  were  obtained,  and  other  diverse 
operational  parameters  involved. 

To  the  authors  knowledge,  no  experiment  has  been 
performed  to  measure  drag  on  a sphere  in  a Poiseuille  flow 
with  a device  having  a sensitivity  of  0.1  micrograms. 
Together  with  this  precise  thermal-microbalance  an  accurate 
gas  flow  meter  with  an  accuracy  of  ±1%  is  employed.  Very 
smooth  teflon  balls  of  different  diameters  (d=  7.88  mm,  9.49 
mm,  11.06  mm  and  12.67  mm)  were  centrally  hung  by  means  of  a 
nylon  string.  The  whole  experiment  was  performed  with  the 
gas  flowing  through  a 22.45  mm  diameter  cylindrical  tube. 


The  drag  force  exerted  on  a sphere  bounded  in  a 
Poiseuille  flow  is  a function  (Figure  4-1)  of  approach 
velocity  Uq,  fluid  density  p,  fluid  viscosity  n,  sphere 
diameter  d,  the  tube  diameter  D,  the  eccentricity  n and  the 
surface  roughness  of  the  ball  e;  viz. 


The  teflon  balls  had  a very  smooth  surface  and  were 
accurately  hung  at  the  center  of  the  cylinder.  This 
eliminated  n and  e and  in  terms  of  dimensionless  parameters 
the  above  relationship  can  be  expressed  as: 


4.3  Coefficient  of  Drag  for  a Bounded  Sphere 


Fo=f{Uo,P,^i,d,D,n,e) 


(4-2) 


2Fq 

PUqA 


(4-3) 
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where  A is  the  cross  sectional  area  of  the  sphere  at  its 
equatorial  plane. 


Figure  4-1  Parameters  effecting  drag  on  a sphere  within  a 

circular  tube 


The  drag  forces  exerted  on  the  sphere  can  also  be 
written  as  : 

2F^=Copu^A  (4-4) 

From  Equations  (4-3)  and  (4-4)  the  drag  coefficient  is  given 
by: 

Ca=f(i?e,p)  (4-5) 

where  Re  is  the  Reynolds  number  based  on  the  mean  approach 
velocity  and  the  tube  diameter  D. 

Klyachko  (1934),  worked  on  sedimentation  of  solid 
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particles  in  fluids  and  gave  the  following  relationship  for 
Cj3  and  Re,  Re  being  based  on  centerline  velocity  and  sphere 
diameter. 


(4-6) 


where  Re  < 1000  and  0.224  £ d/D  £ 0.548  . 

Fayon  and  Happel  (1960)  perfomned  experiments  within  Re 
of  40,  and  Re  had  the  same  definition  as  Klyachko's: 


W ^ 1 + (^-1) 

6n(^ua  i-2.105(  — )+2.087(£)^ 

R R 

where  W is  the  weight  of  the  solid  body,  u is  the  mean 
approach  velocity,  /i  is  the  fluid  viscosity,  a is  the  radius 
of  the  spherical  body,  R the  radius  of  the  enclosing 
cylindrical  pipe,  the  drag  coefficient  for  the  sphere  in 
an  unbounded  system  as  obtained  from  experiment  and  Cg  the 
drag  coefficient  according  to  the  Stoke 's  law. 

In  the  range  10^  < Re  < 10^,  McNown  and  Newlin  (1951) 
studied  the  cylindrical  wall  effect  on  the  drag  of  a sphere 
by  measuring  the  pressure  distribution  of  the  flow  of  air 
around  a sphere.  McNown,  Lee  and  McPherson  (1948)  had 
earlier  related  Cq  with  d/D  as 


(4-8) 


in  the  range  Re  < 200.  They  then  interpolated  the  results  to 
plot  the  variation  of  Cp  with  Re  in  the  range  200  < Re  <10^. 
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Here  too  the  Reynolds  number  (Re)  is  based  on  the  centerline 
velocity  and  the  diameter  of  the  sphere. 

Young  (1960)  obtained  an  expression  for  Cq  by 
suspending  the  spheres  eccentrically  in  a tube.  For  a 
Reynolds  number  between  35  and  1115  he  proposed  a 
relationship 

_ 1 

C[)=Re  ^[94.5(  ^)2*^’^+10]  . (4-9) 

In  the  above  correlation,  Young  defined  Reynolds  number  on 
the  basis  of  the  mean  velocity  of  flow  in  the  pipe  and  the 
pipe  diameter. 

Eichhorn  and  Small  (1964)  investigated  the  fluid 
dyncuaic  forces  on  spheres  suspended  in  Poiseuille  flow. 

Small  spheres  were  suspended  in  a cylindrical  tube 
containing  a laminar  upward  flow  of  water.  Experiment  was 
performed  in  the  Reynolds  number  range  of  80  to  250.  They 
could  not  come  up  with  an  accurate  and  definitive  statement 
about  the  form  of  Re^  and  Cq,  because  of  insufficient  and 
inaccurate  data. 

Clift  et  al.  (1978)  combined  the  relationship  for  many 
observed  data  that  fell  in  the  laminar  range  but  were  mostly 
influenced  by  the  work  of  McNown  and  Newlin  (1951).  As  a 
result,  they  observed  that  for  Re  (based  on  the  approach  or 
centerline  velocity  and  sphere  dicuneter)  lying  between  100 
to  104 
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(4-10) 


where  Cq„  is  the  drag  coefficient  on  a sphere  in  an 
unbounded  flow. 


Figure  4-2  Curvilinear  coordinate  system 

4.4  Fluid  Flow  Analysis 

An  exact  solution  of  the  Navier-Stokes  equation  allows 
the  evaluation  of  fluid  flow  parameters  throughout  the  flow 
domain.  In  our  case,  as  in  most  of  the  practical  instances, 
unfortunately  no  exact  analytical  treatment  of  the  momentum 
equation  is  possible  and  recourse  has  to  be  taken  to  either 
numerical  or  experimental  techniques  to  estimate  the  flow 
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conditions.  In  the  case  of  an  unbounded  sphere  it  is  found 
that  for  a Reynolds  nxiinber  of  3000,  the 

boundary  layer  approximation  to  the  Navier-Stokes  equation 


gives  a correct  estimate  of  the  separation  angle.  On  the 
same  lines  it  is  assumed  that  for  a very  small  dicimeter  of 
the  sphere  compared  to  the  diameter  of  the  cylinder  and  for 
high  flow  rates,  the  boundary  layer  analysis  may  be  able  to 
predict  the  separation  angle. 

The  steady  state  boundary  layer  equations,  according  to 
Boltze  (1908),  are 


Here  a curvilinear  system  of  co-ordinate  is  adapted, 
(Figure  4-2)  where  x is  measured  along  the  body  surface  and 
y perpendicular  to  it.  Variables  u and  v are  the  velocities 
in  the  x and  y direction,  respectively,  p is  the  pressure,  v 
the  kinematic  viscosity  of  the  fluid,  and  r the  radius  of 
curvature  of  the  body.  Equations  4-11  and  4-12  have  been 
solved  for  an  unbounded  fluid  by  Schlichting  (1968).  We  find 
that  the  boundary  conditions  governing  the  boundary  layer 
equation  over  the  sphere  are  the  same,  both  for  bounded  and 
unbounded  fluid.  At  y = 0 , i.e.  at  the  body  surface,  a no- 


/ 


(4-11) 


d(ur) ^d(vr) 

9y“ 


(4-12) 
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Figure  4-3  Accelerated  flow  along  a sphere  in  a 

tube 


slip  condition  exists  viz.  u=v=0.  Aty=5  (boundary 
layer  thickness)  for  an  unbounded  flow  du/dy  = 0 since  there 
is  potential  flow  beyond  the  boundary  layer  thickness.  In 
our  case,  however,  at  y = 6,  u = Ujj,^(x)  and  hence  du/dy  = 0 
(Figure  4-3).  Thus,  both  the  upper  and  lower  limits  for  the 
boundary  layer  are  satisfied  by  the  same  conditions  and  we 
can  carry  out  the  analysis  as  done  for  a sphere  in  infinite 
flow.  The  change  that  needs  to  be  incorporated  is  the 
variation  of  Ujj^^jj(x)  over  the  body  surface.  It  is  not 
accelerated  as  given  by  the  potential  flow,  but  accelerated 
due  to  the  change  in  the  flow  cross  sectional  area  viz. 
[7T.R2(1  - B^sin^<f>)]  as  well  as  due  to  the  body  curvature. 
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Here  R is  the  radius  of  the  cylinder,  13  is  the  diameter 
ratio  of  the  sphere  and  the  cylinder  and  (p  is  the  angle  from 
the  front  stagnation  point  at  which  the  cross  sectional  area 
is  being  measured.  As  a result  will  vary  as 

UoSin<^/( l-B^sin^c^) . 

The  computations  are  carried  out  by  choosing 

• • (4-13) 


andl!{x,y)=  ^ [Uixfi(ti) +2u3X^f3(Ti) +3u5x5f5(il) +. . . (4-14) 

\ 2ui 


where  u^^,  U3,  Ug  etc.  are  the  unknown  coefficients  for 
dependence  on  x;  f2(ii),  f3(Tj),  f5(T|)  etc.  are  the 
functional  coefficients  as  obtained  from  the  Blasius  series; 
T|  = y \/(2u3/v)  is  the  dimensionless  distance  from  the  wall 
with  V being  the  kinematic  viscosity  of  the  fluid  and 
T(x,y)  is  the  stream  function.  For  an  axially  symmetric  case 
we  have  for  the  x-direction  velocity 


u(x) 

} 


^ 1 3(i|ir)  _ 5i|f 
r dy  "^y 


(4-15) 


which  from  Equation  4-14  gives 


U(X)  =UiXf^  +2U2X^f^  +3UgX^f5  +4u-jX^f^  ...  ( 4-16 ) 

If  the  outer  velocity  is  accelerated  according  to  our 
assumption: 
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(4-17) 


where  K is  an  arbitrary  constant,  then 


/ 


(4-18) 


where  the  terms  in  the  parentheses  are  obtained  by  expanding 
(sin<^)/(  1-J3^sin^<^)  as  a power  series  and  the  terms  T3, 

T5,  T7  etc  are  computed  and  found  to  be 


For  a sphere  in  infinite  flow  it  is  observed  that  J3=0  and 
K=1.5. 

Since  (f>={K/r) , where  r is  the  radius  of  curvature  at  a 
distance  x from  the  front  stagnation  point,  we  can  compare 
the  coefficients  of  x in  Equations  4-13  and  4-18  to  get: 

Ui=  K UoT^/r  , U3=  K Uq  T3/r^  , U5  = K Uglg/r®  . . . and  so  on. 
Substituting  the  values  of  these  coefficients  into  Equation 
4-16  we  obtain: 


T3  = 6J32  - 1 

T5  = 120fi^  - 6OJ32  + 1 

T7  = 5040B®  - 420013'^  + 546J32  - 1 


(4-19) 


The  shear  stress  is  obtained  at  the  surface  (y=0)  by 


Newton's  equation: 
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(4-20) 


Since  Tf=yJ  ( 2u^l/v  ) , we  get 


/ du . _ , du  dq 


du  . 

o^ 


2KUq 

rv 


(4-21) 


As  a result  of  this,  the  shear  stress  is  now  given  by 


rvTT  r / X > ^//  . o , X . 3^// 

To=A^^  +2(_)  f3 


To  V c / / 

_+3(  _ )^fc  — 
3!  ' r ^ ^5! 


+•••}] 


(4-22) 

In  order  to  make  a comparison  with  the  unbounded  flow  case, 
we  set  K=1.5  and  obtain 


7.3485 


O.SpU^ax 


(4-23) 


where  Tq  = the  shear  stress 

p = density  of  the  fluid 

Umax  ~ maximxim  velocity  in  the  conduit 

Re  = flow  Reynolds  number  based  on  sphere  diameter 


„ 2<p^T2fi'  3<l>^T^fi^  40'^r'^f7^ 

and,  F(^,cp)  = ^ h ^ 31  ^ ^ (4-24) 

values  of  fi”,  £3”/  £5"  and  f^”  can  be  obtained  from  the 
literature  and  are  provided  in  Appendix  A. 

For  flow  separation  Tq  = 0 i.e.  F(fi,0)=  0.  Substituting 
the  relevant  values  of  and  f'^  (i=l,3,5  ...)  we  get: 
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<^(0.9277)  +<^^(2. 095p2 -0.3641 ) +0^(3. 2466  P'^-1. 5495  p2)  + 
<^'^(4.396P®-3.4907p^+0.4687p2-0.0015)  =0  (4-25) 

where  $ is  the  separation  angle.  We  find  that  the  angle  of 
separation  is  dependant  on  the  sphere  diameter,  if  the 
cylindrical  wall's  dimension  is  unchanged,  unlike  in  the 
case  of  unbounded  fluid  where  the  angle  of  separation  is 
constant.  Table  4-1  provides  values  of  flow  separation  angle 
as  a function  of  sphere-to-cylinder  diameter  ratio  (J3=d/D). 


Table  4-1  (B=d/D)  vs.  separation  angle  {4>sep) 


Diameter  ratio  (13) 

Separation  Angle  (<^gep)  degrees  I 

0.0 

109.60  1 

0.1 

123.86  1 

0.2 

109.36  1 

0.3 

108.46  1 

0.4 

92.61  1 

0.5 

1 

LO 

ro 

• 

00 

00 

0.6 

90.53  1 

0.7 

106.16  1 

Figure  4-4  shows  the  variation  of  shearing  stress  over 
the  surface  of  the  bounded  sphere  for  a laminar  boundary 
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Angle  from  front  stagnation  point  (degrees) 


Figure  4-4  Shear  stress  distribution  on  a bounded  sphere  for 

laminar  boundary  layer 
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Dimensionless  distance  (r|) 


Figure  4-5  Velocity  distribution  as  a function  of  r) 


layer.  The  points  at  which  the  shear  stress  takes  null 
values  are  the  front  stagnation  point  (0=0)  and  the 
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separation  point  (0=0gg_) . Comparison  of  the  sphere  in  free 
flow  with  different  spheres  in  bounded  flow  is  afforded  by 
choosing  K=1.5.  For  B=0,  the  separation  angle  gives  the 
exact  value  as  has  been  provided  by  Schlichting  (1968). 
Figure  4-5  gives  a plot  of  the  velocity  distribution  for 
various  spheres  as  a function  of  the  dimensionless  abscissa 
T|.  The  results  have  been  plotted  at  a point  where  the  flow 
angle  from  the  front  stagnation  point  is  30  degrees.  This 
value  of  angle  is  based  on  the  following  reason.  It  has  been 
found  that  the  boundary  layer  approximation  on  an  unbounded 
sphere  takes  the  potential  flow  as  the  outer  condition  and 
from  experiments  it  is  found  that  this  assmnption  is  true 
until  an  angle  of  30  degrees.  It  is  reiterated  that  boundary 
layer  analysis  gives  a realistic  result  for  the  case  of  an 
unbounded  sphere  only  when  the  Reynolds  number  is  around 
3000.  It  is  assumed  here  that  at  higher  flow  rates  and  for  a 
small  diameter  of  the  sphere  with  respect  to  the  cylinder, 
the  above  analysis  may  shed  some  light  on  the  separation 
angle.  However  in  order  to  get  some  accurate  results  in  this 
area,  numerical  analysis,  should  be  employed.  Chapter  6 is 
devoted  to  the  solution  of  the  Navier-Stokes  equation,  for 
the  geometry  of  interest  to  this  research,  by  finite 


difference  methods. 


CHAPTER  5 

MODELLING  OF  EXPERIMENTAL  DATA 
5 . 1 Error  Analysis 

During  any  experiment  a collection  of  data  is  usually 
obtained  for  different  items.  The  data  acquired  for  each 
individual  parameter  then  combine  to  provide  the  actual 
result.  Thus  in  order  to  determine  the  coefficient  of  drag, 
we  must  have  the  values  for  the  drag  force,  flowrate, 
diameters  of  the  sphere  and  the  cylinder  and  the  value  for 
the  density  of  the  gas.  Each  of  these  parameters  has  some 
estimated  or  characteristic  error  associated  with  it  and 
they  combine  to  give  the  uncertainty  in  the  final  result.  In 
general  a quantity  x which  is  a function  of  u,v,w  ...  has  an 
error  associated  with  it  which  can  be  briefly  stated  as 
follows  (Bevington,  1969): 

X = f(u,v,w  ...)  (5-1) 

The  variance  for  x in  terms  of  variances  a\,  ...  which 

are  actually  known  is 


+ Ov( 


^)  2+2o^ 
OV 


uv 


du  dv 


(5-2) 


The  covariance  can  be  neglected  on  the  following 
grounds.  If  the  fluctuations  in  u,v  are  uncorrelated,  then 
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one  hopes  to  obtain  almost  equal  number  of  positive  and 
negative  values  for  this  term  and  expect  its  contribution  to 
vanish.  On  the  basis  of  the  above  mathematical  background 
one  can  now  proceed  to  perform  error  analysis  for  the 
experimental  values  obtained  in  the  current  study. 

If  'F'  is  the  drag  force  on  the  sphere,  ' p'  the  density 
of  the  gas,  'U'  the  mean  velocity  of  flow,  'A'  the 
characteristic  area  and  Cq  the  coefficient  of  the  drag  then 


2 


=KU^AC^ 


(5-3) 


where  K = p/2 

In  the  experiment  the  value  of  the  flow  rate  was 
measured  and  velocity  was  deduced  from  it.  On  substituting 
the  values  of  mean  velocity  'U'  and  area  'A'  in  terms  of 
flow  rate  and  diameters,  i.e.  U = (4Q)/(7tD^)  and  A = 
(7Td2)/4,  where  'Q'  is  the  flow  rate  and  'D'  and  'd'  are  the 
diameters  of  the  cylinder  and  sphere  respectively  we  get 

F=KQ{^)^Co  (5-4) 


(5-5) 


To  get  the  final  variance  through  Equation  (5-2)  one 
needs  to  get  the  partial  derivatives  from  Equation  (5-5)  and 
they  are 


_1  2D  F 
dD  K Q 


(5-6) 
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5Q3__  2 F 
dd  K Q 

^Q)_  1 1 2)2  1 

dF  d Q 

?Sr  = -1  (2)2  JL 

do  K d q2 


(5-7) 

(5-8) 

(5-9) 


Therefore,  following  Equation  (5-2)  and  neglecting  the 
covariance,  we  get 


(5-10) 


a (221)2^^  i2l£)2^2z  (211)2^^  {2lJL)2 

Co  ^2  ^2  q’  k2  ^3  k2  q’  j^2  j2  q2 


(5-11) 


Dividing  Equation  (5-11)  throughout  by  C^,  we  get 


OcJ  _ . On. . . Ow. . . Oc.. . a 


r-2 


= 4(-^)2+4  (_I^)2+(^)2+(^)2 

D d F C) 


(5-12) 


In  the  current  experiment  the  drag  force  'F'  can  be 
measured  with  an  uncertainty  of  ±0.001%,  the  flowrate  'Q' 
has  an  uncertainty  of  ±1%  while  the  uncertainties  in  the 
sphere  diameter  and  cylinder  diameter  are  ±0.1%  and  ±2.0% 
respectively.  These  are  the  maximum  errors  prescribed  by  the 
manufacturers.  As  a result,  the  maximum  variance  we  get  is 
(acD/Co)=  ±0.041279. 

If  eccentricity  of  the  sphere  is  also  considered  as  one 
of  the  governing  factor,  f(13),  then  it  appears  as  f(B)2  in 
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the  error  propagation  equation.  If  the  error  due  to 
eccentricity  is  1%  then  (Oj^q/Cq)  = ±0.0439137.  On  inspection 
it  is  evident  that  the  major  error  in  this  experiment  is  due 
to  the  measurement  of  cylinder  diameter.  It  is  very  easy  to 
obtain  precisely  bored  tubes.  Moreover  as  far  as  error  in 
the  flow  is  concerned,  it  should  be  noted  that  the  range  of 
the  instrument  was  0-50  Lt/min.  The  maximum  error  occurs  at 
the  higher  ranges  of  flow.  While  conducting  viscosity 
experiments,  one  has  to  employ  the  lower  flow  rates  as  has 
been  discussed  ahead.  As  such  while  designing  a 
sophisticated  drag  based  viscometer,  some  of  the  typical 
error  values  that  can  be  safely  assumed  for  now  are: 
±0.0001%(for  drag  force),  ±0.1%  (for  flow  rate),  ±0.001% 

(for  sphere  diameter),  ±0.01%  (for  cylinder  diameter)  and 
±0.01%  (for  eccentricity).  These  combine  to  provide  a total 
uncertainty  in  drag  coefficient  measurement  as 

( ^cd/^d ) . 003 . 

5.2  Data  Modelling 

An  efficient  way  to  model  the  drag  force  data,  is  to 

plot  the  variation  of  coefficient  of  drag  (Cq)  with  the  flow 

Reynolds  number  (Re).  Since  the  flow  characteristics  vary 

with  the  Reynolds  number,  a convenient  way  is  to  model  the 

data  for  each  flow  phenomenon.  Before  one  starts  on  this 

venture,  one  should  be  careful  in  ascertaining  the  Reynolds 

number  which  by  definition  is  the  ratio  of  the  fluid's 

« 

inertial  force  to  the  viscous  force  it  exerts.  The  inertial 
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force  is  a function  of  the  flow  velocity.  For  an  infinite 
flow  over  a sphere  this  velocity  is  the  constant  approach 
velocity.  But  for  a flow  in  a pipe  when  the  fully  developed 
condition  has  been  reached,  the  velocity  distribution  is 
parabolic  and  is  given  by  the  Hagen-Poiseuille  formulation 

(5-13) 

Therefore  when  comparing  the  bounded  flow  with  an  infinite 
flow  to  discern  the  wall  effects,  the  disparity  in  the 
velocity  must  be  appropriately  taken  into  account. 


Figure  5-1  Velocity  profile  seen  by  different  spheres 


• • 
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Also,  as  can  be  seen  from  Figure  5-1, the  two  different  sized 
spheres  encounter  different  velocity  distributions  for  the 
same  cylinder.  While  performing  dimensional  analysis  authors 
have  generally  chosen  two  different  characteristic 
velocities  to  calculate  the  Reynolds  nximber.  One  is  the  mean 
velocity  in  the  tube  [Uj„ggjj=Q/ (ttR^  ) where  Q is  the 

flowrate  and  'R'  is  the  tube  radius]  while  the  other  is  the 
centerline  velocity  • The  choice  of  is 

based  on  the  grounds  that  the  velocity  distribution 
preceding  the  particle,  which  is  a function  of  pipe  Reynolds 
number,  will  have  considerably  more  influence  on  the  drag 
coefficient.  But  the  use  of  looks  erroneous  on  the 

grounds  that  not  all  spheres  face  the  mean  velocity.  The 
smaller  sphere  in  Figure  5-1,  encounters  more  of  the 
velocity  near  the  centerline  which  is  much  more  than  the 
mean  velocity.  Figures  5-2  and  5-3  show  the  variation  of  Cq 
with  Re  and  sphere  to  cylinder  size  (d/D)  respectively,  as 
has  been  obtained  in  the  current  experiment.  By  choosing 
Umean  define  the  Reynolds  number,  a plot  in  the  form  of 
Figure  5-2  is  obtained.  Figure  5-3  shows  that  for  B( sphere 
to  cylinder  diameter  ratio) =0.5  the  coefficient  of  drag  (Cq) 
is  minimum.  This  can  be  explained  on  the  basis  of  the 
velocity  distribution  encountered  by  the  spheres.  The 
smaller  sphere  does  not  lie  in  the  path  of  the  mean  velocity 
and  faces  more  of  the  velocity  grouping  around  the 
centerline.  This  velocity  distribution  is  larger  than  the 
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V 

REYNOLDS  NUMBER 


Figure  5-2  Coefficient  of  drag  (Cq)  vs.  Reynolds  number  (Re) 

for  different  spheres  in  a tube 
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Figure  5-3  Coefficient  of  drag  (Cd)  vs.  diameter  ratio  (d/D) 

for  different  Reynolds  number 
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mean  velocity.  Thus  choosing  as  the  characteristic 

velocity  for  the  expression 


C 


^ DRAG 


(5-14) 


overestimates  the  coefficient  of  drag,  since  the  actual  drag 
force  exerted  is  due  to  a larger  velocity  than  the  mean 
velocity.  The  reverse  is  true  in  the  case  of  larger  spheres. 
The  use  of  Uj„ax  more  justified  since  all  the  spheres 
which  are  centrally  suspended,  at  least  meet  it.  Figure  3-3 
had  made  a comparison  of  the  coefficient  of  drag  as  obtained 
by  various  authors  on  the  basis  of  centerline  velocity. 
However,  a more  appropriate  choice  of  velocity  can  be  had. 
This  is  based  on  the  ground  that  there  is  a velocity 
distribution  flowing  over  the  sphere  and  a characteristic 
velocity  must  be  derived  from  it. 

When  obtaining  the  coefficient  of  drag  (Cq)  from 
Equation  4-1,  one  observes  that  it  is  the  ratio  of  the  drag 
force  (D)  experienced  by  the  body  to  the  dynamic  force 
carried  by  the  fluid  (0.5(17^).  So  a more  judicious  choice  of 
the  characteristic  velocity  would  be  one  that  exerts  the 
same  dynamic  force  as  enthrusted  by  the  velocity 
distribution  seen  by  the  sphere.  The  velocity  distribution 
in  the  pipe  is  given  by  the  Hagen-Poiseuille  Equation  5-13. 
The  total  dynamic  force  exerted  by  this  velocity 
distribution  on  the  sphere  is  (1/2)  where 
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the  integral  is  bounded  between  0 and  'a'  with  'a'  and  R 
being  the  radius  of  the  sphere  and  cylinder,  respectively. 
The  characteristic  velocity  (Uq)  is  given  by 


(5-15) 


Uq  gives  the  same  dynamic  force  as  exerted  by  the  velocity 
distribution  seen  by  the  sphere.  Since  the  velocity 
distribution  is  the  parabolic  profile  given  by  Equation  5-13 
we  get  from  Equation  5-15 


Un=2U 


mean 


(5-16) 


Hence,  the  characteristic  velocity  (U^)  is  equal  to 

V^(  l-2fi^/3+6^/5 ) , where  J3=r/R  and  is  the  centerline 

velocity  and  twice  With  this  characteristic  velocity 

it  is  fruitful  to  perform  an  analysis  between  a bounded  and 
an  unbounded  flow.  In  the  derivation  of  the  characteristic 
velocity,  the  average  value  obtained  was  for  a linear 
distribution.  It  is  obvious  that  a better  averaging  for  this 
type  of  flow  would  be  by  taking  the  cross-sectional  average. 
If  we  denote  the  characteristic  velocity  obtained  in  this 
way  by  then 
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a 


u^lnrdr 


(5-17) 


2rrrdr 

D 


and  it  provides 


(5-18) 


For  the  same  flow  rate,  the  velocity  has  a constant  value  Uo, 
for  infinite  flow  and  or  for  a bounded  flow.  The 
result  is  that  the  Reynolds  number  in  an  infinite  flow  and 
bounded  flow  is  now  defined  by  choosing  the  characteristic 
velocity  and  is  given  by  Re=(  pU^d)//Li  where  d is  the  sphere 
diameter,  and  is  given  by  the  term  l-2B^/3+J3'^/5 ) 

where  is  the  centerline  velocity  and  B is  the  diameter 

ratio  between  the  sphere  and  cylinder.  Equipped  with  these 
definitions  for  the  Reynolds  number  and  the  characteristic 
velocity,  a comparison  between  bounded  and  unbounded  flow 
brings  out  remarkable  closeness  in  the  values  of  the 
coefficient  of  drag  for  different  cases.  The  equation  which 
relates  the  value  of  the  coefficiemnt  of  drag  for  bounded 
spheres  (C^q)  as  is  observed  in  this  work,  to  the  published 
values  of  coefficient  of  drag  for  unbounded  spheres  ( Cqoo)  is 
given  by 


(5-19) 


where  AfP  when  we  choose  Uq  and  A?=0.9pP  when  is  chosen. 
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The  error  is  given  by  : ERROR  = [ (Cd„  - Cdj3)/Ci3„|x  100  where 
Cqq,  can  be  obtained  from  standard  relationships  as  given  by 
Equations  2-3  and  2-4.  Table  5-1  provides  the  absolute 
values  for  error  in  drag  coefficient  evaluation.  The  error 
is  obtained  when  comparing  our  experimental  data  to  the 
values  obtained  for  Cqq  from  Equation  5-19.  In  Table  5-1  El 
and  E2  are  respectively  the  absolute  error  (%)  obtained  by 
choosing  Uq  and  It  should  be  noted  that  although  the 

maximiim  error  is  around  6.5%  the  relationship  for  C^a,  has 
itself  an  error  of  more  than  6.5%  associated  with  it.  This 
has  been  noted  by  Clift  et  al.  (1978)  who  obtained  the 
relationship  (Cpo,)  for  the  unbounded  sphere  by  fitting 
experimental  data  through  linear  regression.  The 
experimental  data,  as  collected  during  this  study,  was  first 
curve  fitted  through  least  square  method  by  assuming  a 
polynomial  form  of  the  form: 

y=  a + bx  + cx^  + dx^  + ex^  + ...  (5-20) 

Table  5-2  shows  the  coefficients  obtained  for  fitting 
experimental  values  of  this  work  for  different  spheres  and 
Reynolds  niimber  between  200  and  800. 


Table  5-1  Table  showing  absolute  errors  in  coefficient 

of  drag  measurement  using  relationship  (5-19) 
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Table  5-2  Coefficients  used  for  fitting  experimental  data 


J3 

a 

b 

c 

d 

r 

0.35 

1.3651 

-3.662 

E-3 

6.4061 

E-6 

-4.265 

E-9 

0.99 

0.98 

0.42 

1.8944 

-6.052 

E-3 

9.9561 

E-6 

-5.710 

E-9 

0.99 

0.99 

0.49 

1.5271 

-3.418 

E-3 

4.5935 

E-6 

-2.232 

E-9 

0.99 

0.99 

0.56 

2.2202 

-5.979 

E-3 

8.1808 

E-6 

-3.773 

E-9 

0.99 

0.98 

The  coefficients  thus  obtained  allowed  us  to  perform 
comparison  at  different  Reynolds  number.  Although  this  form 
of  comparison  can  be  extended  beyond  Reynolds  number  650  but 
it  has  not  been  reported  here  due  to  lack  of  data  at  higher 
Reynolds  number  for  B=0.35  and  0.42.  Also,  below  Reynolds 
number  of  250,  the  data  were  tested  and  the  results  were 
consistent  until  Reynolds  number  of  200.  Below  the  Reynolds 
number  of  200,  the  fluid  flow  phenomenon  is  different,  and 
data  from  this  work  was  not  sufficient  enough  to  provide  a 
new  correlationship. 


CHAPTER  6 

NUMERICAL  STUDY  of  FLUID  FLOW 
6 . 1 Introduction 

There  are  several  tools,  in  the  foinm  of  niomerical 
algorithms,  available  at  present  to  perform  computational 
studies  of  the  fluid  momentum  equation.  According  to  the 
complexity  of  the  physical  phenomenon  the  Navier-Stokes 
equation  can  either  be  simplified  to  a boundary-layer  type 
equation,  parabolized  Navier-Stokes  equation,  or  can  be 
retained  in  its  full  form.  A comprehensive  discussion  on  the 
solution  of  these  fluid  mechanics  equations  by  finite 
difference  method  has  been  provided  by  Anderson  et  al . 

(1984).  The  full  blown  Navier-Stokes  equation  itself  can  be 
differently  treated  for  compressible  or  incompressible  fluid 
flow.  Numerical  techniques  that  handle  compressible 
situations  run  into  convergence  difficulties  for 
incompressible  case.  As  a result,  the  treatment  of  these  two 
flow  situations  require  separate  methodologies.  For  the 
circumstances  that  govern  our  problem,  viz.  laminar  flow  of 
gas  over  a cylindrically  surrounded  rigid  sphere,  the 
incompressible  Navier-Stokes  equation  has  been  solved.  There 
are  two  popular  approaches  to  handle  the  incompressible 
fluid  momentxim  equation:  a)  Stream  Function-Vorticity  method 
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and  b)  Primitive-  Variable  approach.  The  application  of  the 
former  method  has  been  successful  for  two  dimensional  cases 
but  runs  into  severe  difficulties  in  three  dimensional 
problems  and  during  pressure  calculations.  In  spite  of  this 
the  Streamfunction-Vorticity  method  has  some  distict 
features  to  offer.  On  the  other  hand  Primitive-  Variable 
approach  is  widely  used  these  days  and  is  not  limited  to  two 
dimensional  situations.  In  the  current  study,  both  the 
methods  have  been  employed. 


The  Navier-Stokes  equation  in  a cylindrical  coordinate 
(r,0,z)  system  has  been  used  for  the  computational  fluid 
dynamics  study  (Figure  6-1). 


6.2  Governing  Equations  and  Discretization 


Y 


Figure  6-1  Cylindrical  coordinate  system 


74 


The  flow  is  assumed  to  have  symmetry  about  the  '0'  direction 
and  as  a result  we  are  interested  in  only  axisymmetric  flow. 
Strictly  speaking  this  narrows  the  study  to  Reynolds  numbers 
where  asymmetry  has  not  set  in.  For  a sphere  in  infinite 
flow,  this  symmetric  flow  restriction  limits  the  Reynolds 
number  to  values  within  240. 

6.2.1  Primitive  Variable  Equations 

If  we  assume  a Newtonian  fluid  with  constant  viscosity 
and  the  absence  of  body  forces,  then  for  two-dimensional 
axisymmetric  flow  (i.e.  symmetry  and  no  velocity  component 
in  the  '0'  direction)  the  expressions  for  the  momentum  and 
continuity  equations  are: 


=-  1 dp  du. 


(6-1) 


dv  dv  dv  __  1 dp  . d^v  d^v  .ldv_  v 

WE  W^  p'J?  ~d^  r'^  72 


(6-2) 


du^l  d(rv) 

W^  7 dr 


(6-3) 


Here  't'  is  the  time,  'u'  is  the  velocity  component  in  the 
'z'  direction,  'v'  is  the  velocity  component  in  the  'r' 
direction,  'v'  is  the  kinematic  viscosity  of  the  fluid,  ' p' 
is  the  density  of  the  fluid  and  'p'  is  the  pressure. 
Equations  (6-1)  and  (6-2)  are  the  momentum  equations  in  the 
'z'  and  'r'  directions  respectively  while  Equation  (6-3)  is 
the  fluid  continuity  equation.  The  above  set  of  equations 
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are  referred  as  the  "Primitive”  equations  since  they  are 
written  in  terms  of  the  primitive  variable  'u','v'  and  'p'. 
6.2.2  Stream  Function-Vorticity  Equations 

The  pressure  term  in  Equations  6-1  and  6-2  can  be 
eliminated  by  cross  differentiation  and  applying  Equation 
6-3  and  the  vorticity  (C)  relation  for  axisymmetric  flow 

the  following  equation  with  vorticity  as  the  dependent 
variable  is  obtained. 


« -z? -V  ( i!£  .i!i  .1 « - 

WE  Wz  Wr  r rWr  j-2 


(6-5) 


If  a slightly  modified  version  of  the  continuity  equation 
(Equation  6-3), 


Pr  1 3(rv)  § d(rv) 

T~WF~ 


(6-6) 


is  added  to  the  left  side  of  Equation  (6-5),  we  get  the 


"vorticity  transport"  equation: 


Wi  dz  dr  0^2  0^2  ^2' 


(6-7) 


Introducing  the  incompressible  stream-function  ( i|r)  defined 
by 
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, l4^  = -v  (6-8) 

r or  r dz 

into  the  continuity  equation  (Equation  6-3)  gives  the 
following  equation  for  the  vorticity  in  terms  of  the  stream- 
function 


Equations  6-7,  6-8,  and  6-9  constitute  the  'Streamf unction  - 
Vorticity'  equations. 

6.2.3  Discretization 

In  order  to  solve  the  partial  differential  equations 
numerically,  first  the  physical  space  is  divided  into  a 
discrete  number  of  points  which  conform  to  some  cartesian 
system  and  form  a grid.  Then  the  governing  equations  - 
Equations  6-1,  6-2  and  6-3  for  the  primitive  method  and 
Equations  6-7,  6-8  and  6-9  for  the  stream  function-vorticity 
formulation  are  discretized  to  break  them  down  into 
algebraic  equations  that  can  be  solved  niimerically . Several 
techniques  like  the  use  of  Taylor  series,  polynomial 
fitting,  integral  method,  and  the  control  volume  approach 
can  be  utilized  for  discretization.  The  efficacy  of  one 
method  over  the  other  is  problem  dependant  and  cannot  be 
ascertained  beforehand.  However  the  control  voliime  approach 
offers  some  interesting  features.  It  provides  a balance  of 
some  physical  quantity  on  a region  in  the  vicinity  of  the 
grid  point.  This  means  that  the  governing  physical  law  is 
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satisfied  over  a finite  region  rather  than  a single  point. 
Moreover,  the  different  equations  formed  by  the  control 
volume  approach  tend  to  have  the  conservative  property  which 
is  essential  to  many  problems. 

An  inspection  of  the  transport  equations  shows  that  the 
dependent  variables  (u,v,C)  obey  a certain  conservative 
principle.  In  general  if  is  a dependent  variable,  then 

the  differential  equation  in  a conservative  form  can  be 
written  as 


dz  9r  ^ ' 


(6-10) 


where  'T'  is  a diffusion  coefficient  and  'S'  is  the  source 
term.  In  Cartesian  tensor  notation  we  have 


d(p<p)  ^ d 30 

dxl 


(6-11) 


where  is  the  fluid  velocity  in  the  direction  The  four 
terms  appearing  in  the  equation  are  the  unsteady  term,  the 
convective  term,  the  diffusive  term,  and  the  source  term. 

The  equation  as  a whole  represents  the  conservation  of  flux. 
The  total  flux  consists  of  a)the  convective  flux  (fU^$) 
due  to  the  motion  of  the  fluid  and  b)the  diffusive  flux 
[r(3$/3x^)]  generated  by  the  gradients  of  $.  The  physical 
meaning  of  expressing  the  terms  in  such  an  order  has  been 
given  by  Patankar  (1980).  One  can  easily  see  that  Equations 
6-1, 6-2  and  6-7  can  easily  be  cast  into  the  generalized  form 
as  given  by  Equation  6-10.  Therefore  it  is  meaningful  to 


78 


discretize  the  generalized  equation  and  then  substitute  and 
solve  for  the  variables  of  interest  like  u,v  or  5* 


In 


Figure  6-2  Control  Volume  around  a grid  point  'P' 

T- 

In  the  control  volume  approach  the  discretization  is  done  by 
integrating  the  governing  equation  over  a region.  In  Figure 
6-2,  a control  volume  is  chosen  around  the  grid  point  'P' 
and  Equation  6-10  is  integrated  over  it.  Discretization  in 
this  manner  results  in  the  following  algebraic  equation  (see 
Appendix  A and  Patankar  (1980)  ): 

ap$p=ag$E'*'^w^w+^N®N^®s^s"*’^  ( 6—12 ) 

The  definition  of  the  terms  ap, a^, a^^, a^,, ag  and  b is  given  in 
Appendix  A.  This  resulting  algebraic  equation  is  solved 
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numerically  by  observing  boundary  conditions,  the  nature  of 
the  equation,  and  applying  a proper  differencing  scheme  on  a 
proper  set  of  grid  lines.  The  next  task,  therefore,  is  to 
generate  an  appropriate  grid. 

6 . 3 Grid  Generation 

The  generation  of  grid  points  within  the  physical  space 
is  an  important  consideration.  More  often  than  not,  results 
are  needed  on  surfaces  of  the  problem  and  we  must  have  grid 
points  on  it.  From  Figure  6-3  we  observe  that  if  grid  lines 
are  aligned  along  the  z and  r direction  of  the  cylinder  then 
it  is  difficult  to  ascertain  values  on  the  sphere  surface. 
For  example,  points  A,C  do  not  lie  at  the  intersection  of 
grid  lines  like  points  B and  D.  Therefore,  some 
interpolation  needs  to  be  done  to  calculate  values  at  points 
like  A and  C.  On  the  other  hand,  if  we  choose  a spherical 
coordinate  system  then  the  cylinder  surface  has  too  few  grid 
points  on  it  and  the  cylinder  is  not  accounted  for  properly. 
This  poses  a classic  problem  where  one  orthogonal  coordinate 
system  (spherical)  is  embedded  in  another  (cylindrical).  If 
the  size  of  the  sphere  is  very  small  compared  to  the  tube 
diameter,  then  a)  a cylindrical  system  would  have  given  some 
fruitful  results  provided  no  values  were  needed  on  the 
sphere  surface  and  b)  a spherical  system  would  have  provided 
results  for  the  sphere  but  very  little  and  debatable  results 
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Figure  6-3 . A r-z  coordinate  for  a sphere  in  a cylinder 

for  the  whole  flow  situation.  But  in  the  current  problem, 
the  size  of  the  sphere  is  not  necessarily  small.  Also,  the 
imposition  of  one  orthogonal  cartesian  coordinate  on  the 
physical  domain  will  reguire  some  sort  of  interpolation  for 
the  implementation  of  the  boundary  conditions.  Since 
boundary  conditions  play  a dominant  role,  the  interpolation 
creates  inaccuracies  at  sensitive  places.  In  addition, 
uneven  grid  spacing  near  the  boundaries  are  created  which 
make  the  programming  cvimbersome. 

The  dilemma  may  be  overcome  by  choosing  a body  fitted 
coordinate  system.  As  the  name  indicates,  the  coordinate 
system  is  aligned  along  the  body  surfaces  and  need  not  be 


81 


orthogonal.  This  is  implemented  by  transformation  of  the 
physical  problem  (from  an  unevenly  spaced  and  non-orthogonal 
gridded  region)  to  a computational  space  (evenly  spaced  and 
orthogonal  grids).  A typical  evolution  of  the  computational 
space  from  a physical  domain  is  shown  in  Figure  6-4.  It  is 
to  be  noted  that  the  computational  region  is  obtained  by 
twisting  and  deforming  the  physical  domain.  The  result  of 
this  systematic  deformation  leaves  a space  where  we  have 
orthogonal  grids  at  equal  spacing.  The  fluid  dynamics 
equation  are  now  solved  on  the  computational  domain.  The 
central  issue  in  this  process  is  identifying  the  location  of 
the  grid  points  in  the  physical  domain.  This  means  obtaining 
the  X and  y values  of  the  grid  points  in  the  physical  space 
that  correspond  to  P(i,j)  in  the  computational  domain. 

During  this  identifying  process  for  grid  points  several 
desirable  features  must  be  aimed  at  (Hoffman  1989) 

a)  the  mapping  should  generate  one-to-one  correspondence, 
i.e.  grid  lines  of  the  same  family  should  not  cross 

b)  the  grid  distribution  should  be  smooth 

c)  the  grid  lines  should  be  orthogonal  and 

d)  there  should  be  an  option  for  grid  clustering. 

Some  of  the  promising  techniques  that  help  to  obtain 
some  and  sometimes  all  of  the  above  desirable  features  are: 
the  algebraic  method,  the  partial  differential  method  and 
conformal  mapping.  A detailed  discussion  of  these  can  be 
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elicited  from  Thompson  et  al.  (1985).  For  our  problem  we 
have  used  two  options  for  generating  the  grid  points: 


b.  Computational  domain 


Figure  6-4.  Transformation  from  physical  to  computational 

domain 
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an  axisymmetric  (r-z)  coordinate  and  a body  fitted 
coordinate  with  the  grid  generated  by  the  partial 
differential  method.  Amongst  the  various  options  in  this 
class  (elliptic,  parabolic  and  hyperbolic),  an  elliptic  grid 
generator  was  chosen  since  it  is  widely  used  in  cases  having 
a well  defined  geometry.  If  (x,y)  are  the  variables  of  grid 
points  in  the  physical  space  and  (5/"n)  '^he  computational 
space,  then  the  values  of  (x,y)  within  the  physical  space 
are  determined  by  solving  the  set  of  Equations  6-13  and  6- 
14. 


To  solve  Equations  6-13  and  6-14  by  iterative  methods  like 
the  Gauss-Seidel  technique,  they  are  arranged  as: 


ax^  g -2bX]^^  ■'■^^1111 


(6-13) 


ay55-2b75^+cy^i^=0 

a = + y2 

b = X5X^+  y^y^ 


(6-14) 


where 


c = x|  + y 


b 


(2A^Ati) 


(6-15) 
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(A5) 


(yi+l,j'^y ■^TTTT^  ^Y±,j+l'^y 2 At]  ) 


(Ati) 


( 2 ( Ati  ) 


)} 


(6-16) 

To  begin  the  solution  by  the  elliptic  method,  first  an 
initial  distribution  of  x and  y coordinates  of  grid  points 
within  the  physical  space  is  specified.  The  coefficients  a,b 
and  c are  determined  by  finite  difference  approximation.  A 
second  order  accurate  differencing  scheme  is  employed.  For 
example,  the  value  of  x^  at  interior  points  is  given  by 


’I  2^^i 

and  at  the  boundary  (e.g.,  at  j=l)  it  is  given  by 

^ = -3Xi,i+4Xi,2-Xi,3 
’I  ZSri 

For  the  first  iteration  the  value  of  x and  y are  provided  by 
the  initial  distribution  and  subsequently  from  the  previous 
iteration,  i.e.,  the  computation  of  the  coefficients  lags 
behind  by  one  iteration  level.  This  iterative  process  is 
continued  until  a specified  value  of  convergence  is  reached. 
For  the  current  problem  the  following  strategy  is  adopted 

1=IM1 ,j=JMl 

Y:  abs(xJ;x,.) 

i=2,j=2 


ErrorX = 
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i=IMl ,j=JMl 

ErrorY=  £ 

i=2,j=2 

ErrorT  = ErrorX  + ErrorY 

where  k represents  the  iteration  level . The  convergence  is 
reached  once  ErrorT  £ e,  where  e is  a small  value  and  is 
equal  to  10"^  in  this  case.  On  using  the  r-z  coordinate  the 
sphere  is  approximated  by  a step  like  structure  while  the 
niimerical  grid  generation  exactly  traces  all  the  boundary. 
Figure  6-5  shows  the  grid  lines  in  the  two  situations.  For 
the  current  problem  each  yields  a solution  which  is 
partially  correct. 

6 . 4 Solution  Technique 

Once  the  differential  Equation  6-11  is  discretized  and 
reduced  to  the  algebraic  form.  Equation  6-12  needs  to  be 
solved  to  yield  the  desired  results.  Several  methods  that 
are  employed  to  solve  algebraic  equations  can  be  utilized. 
One  should  note  that  the  resulting  equation  is  two- 
dimensional  and  the  dependence  of  the  coefficients  (ag,ayj 
etc.)  on  the  variable  $ delineates  it  as  a nonlinear 
equation.  Several  iterative  techniques  are  available  to 
solve  this  type  of  nonlinear  equation-  Gauss-Seidel  method, 
Point  Successive  relaxation.  Line  Successive  relaxation, 
etc.  For  our  problem.  Point  Successive  relaxation 
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a.  Physical  region  using  r-z  coordinate 


b.  Physical  region  using  boundary  fitted  coordinate 


Figure  6-5.  Grid  lines  using  different  coordinates 
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was  used  for  the  Stream  function-Vorticity  method  and  Line 
Successive  relaxation  method  for  the  Primitive-Variable 
method.  The  relaxation  parameter  is  introduced  as  follows: 

^ p~^  ^ i'P (6-17) 

where  i=E,W,N  and  S. 

<t>P=^Y.^±<t>±^b  (6-18) 

3p 

If  $p  is  the  known  value  at  'P'  from  previous  iteration  then 
adding  ($p-$p)  to  right  hand  side  of  Equation  6-18  and 
rearranging  terms  one  gets 

<l>P=<l>p-^^('^^i4>i*b-ap<pp)  (6-19) 

3p 

As  the  solution  proceeds,  $p  must  approach  $p.  A relaxation 
parameter  (<o)  is  introduced  to  accelerate  this  process.  The 
value  of  ti)  determines  its  nomenclature;  if  it  is  less 
than  unity  then  it  is  called  underrelaxation,  else  it  is 
overrelaxation . 

The  above  iterative  technique  as  applied  to  Point  and  Line 
Successive  relaxation  can  be  found  in  several  standard  books 
of  Computational  Fluid  Dynamics. 

Before  one  applies  the  iterative  strategy  for  the 
solution  of  the  algebraic  equation,  a generalized  form  is 
formed  so  that  different  differencing  schemes  can  be 
handled.  It  has  been  found  that  for  a convective-diffusive 
equation  such  as  Equation  6-11,  it  is  the  handling  of  the 
non-linear  convective  term  which  assures  the  accuracy  of  the 
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solution.  This  handling  depends  on  the  problem  and  flow 
conditions  and  good  differencing  schemes  try  to  model  them 
accurately.  The  differencing  schemes  can  be  divided  into  two 
groups.  One  is  based  on  one-dimensional  flux  balance  (lower 
order  scheme)  and  cannot  work  well  in  regions  were  flow  is 
not  aligned  with  the  grid  lines.  Some  of  the  well  known 
differencing  methods  falling  under  this  bracket  are  the 
upwind  scheme,  hybrid  scheme,  central  differencing, 
exponential,  power-law,  etc.  The  other  schemes  use  a better 
differencing  technigue  to  handle  grid  to  flow  skewness  and 
some  of  them  are  the  Locally  Analytical  Differencing  Scheme, 


Figure  6-6  Flux  variation  across  a cell  in  Linear  flux 

spline  scheme 
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the  Linear  Flux-Spline  Scheme,  the  Cubic  Flux-Spline  Scheme, 
Controlled  Numerical  Diffusion  with  Internal  Feedback,  etc. 
An  overview  of  these  methods  can  be  found  in  Karki 
et  al.  (1988).  The  problem  with  lower  order  schemes  is  that 
they  do  not  respond  well  to  the  presence  of  sources  or  sinks 
in  the  control  volume  nor  to  multidimensional  effects.  The 
skewness  between  the  grid  and  the  flow  gives  rise  to 
significant  numerical  diffusion  (numerical  error)  and  lower 
order  schemes  fall  short  in  their  performance.  The  higher 
order  scheme  tries  to  remove  these  shortcomings.  The  linear 
flux-spline  scheme  assumes  a linear  variation  of  the  total 
flux,  as  shown  in  Figure  (6-6),  as  opposed  to  the 
preservation  of  flux  that  is  assumed  by  lower  order 
techniques.  In  order  to  combine  all  the  differencing  schemes 
under  one  lunbrella.  Equation  (6-12)  is  rewritten  here  and 
the  coefficients  are  defined  as  follows: 


Sp(pp  -a  g(p£  +a  w<I>w*^n^n'*'^s^s 


where 


ajj  — D jjA  ( 
as  = DgA( 


w 


n 


s 


) + max ( -Fg , 0 ) 
) + max(F„,0) 

) + max(-Fn,0) 

) + max(F_,0) 


The  values  of  and  Fj^  (i=e,w,n  & s)  are  given  in  Appendix 


B. 
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Table  6-1  Value  of  A(jPj)  for  various  schemes 


SCHEME 

FORMULA  FOR  A ( [ P | ) 

CENTRAL  DIFFERENCE 

1 

o 

• 

UPWIND 

1 

HYBRID 

MAX(0,1  - O.SjPj ) 

POWER  LAW 

MAX[0, {1  - 0.5 |P| ) }^] 

EXPONENTIAL 

1P!/[EXP( IPI )-l] 

Physically  is  the  conductance  through  face  i,  is  the 
flow  rate  through  face  i and  P^  is  the  Peclet  number  given 
by  the  ratio  of  the  flow  rate  to  conductance  through  face  i 
i.e.  The  function  A(jPj)  for  lower  order  schemes 

is  given  in  Table  6-1  and  can  be  selected  according  to  the 
method  used.  Unfortunately  for  the  linear  flux-spline  the 
coefficients  cannot  be  tabulated  in  the  table  but  the  values 
of  the  coefficients  are  provided  in  the  works  of  Karki  et 
al.  (1988). 


6 . 5 Implementation 

The  solution  strategy  for  the  Streeim  Function  - 
Vorticity  method  and  Primitive  variable  method  are 
different.  Before  one  starts  working  out  the  details  for 
them  it  is  necessary  to  lay  down  the  boundary  conditions 
(Figure  6-7)  that  regulate  the  solution. 


INFLOW 
(1) 


T.  (2) 


TOP  BOUNDARY  (3) 
(SOLID) 


SOLID  BOUNDARY 


BOTTOM  BOUNDARY  (3) 
(SOLID) 


OUTFLOW 
(4) 


Cl  (2) 


Figure  6-7  Boundary  Conditions  existing  for  the  problem 


Inflow.  For  the  streamf unction  - vorticity  method, two 
different  types  of  velocity  inflow  conditions  were  used 
(surface  1): 

a)  a constant  velocity  at  the  inlet 


r .e. 


u=U 


mean' 


V=0 


/ 


which  yields  the  stream  function  4 = (Ujj^g^j^r^  )/2  and 

b)  a fully  developed  velocity  profile 
i.e.  u=2U^g^j^(  l-r^/R^ ) , which  is  the  Hagen-Poiseuille 
equation.  The  stream  function  becomes 

=2U^gg„[(r2/2)  - (r4)/(4R4)] 
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Both  conditions  a)  and  b)  can  be  put  in  a single  form  like: 

(Ar‘^/4)  + (Br^/2)  (6-22) 

For  condition  a)  A=0  and  and  for  b)  and 

B=Ujj,g^j^.  The  vorticity  in  either  case  is  obtained  from 
Equation  (6-4). 

For  the  primitive  variable  method,  a constant  inflow 
velocity  (u=constant  and  v=0)  was  used. 

Centerline.  The  centerline  is  shown  as  surface  2.  For 
the  streamfunction  - vorticity  method,  the  streamline  can  be 
assigned  an  arbitrary  constant  value,  and  is  taken  as  zero. 
Since  the  streamline  was  given  a constant  value  at  this 
surface,  the  radial  velocity  v becomes  zero  (5T/3z  = 0). 

The  assumption  of  axial  symmetry  results  in  dn/dr  = 0.  With 
dv/dz  = 0 (since  v=0  along  centerline)  and  d\x/dr  = 0 the 
vorticity  is  zero  along  the  centerline. 

In  the  primitive  variable  approach,  a symmetric  boundary 
resulted  in  using  v=0  and  du/dr  = 0. 

Solid  boundary.  The  solid  surfaces  are  the  cylinder  and 
sphere  walls  and  have  been  labelled  as  (3).  The  value  of  the 
stream  function  on  the  walls  of  the  cylinder  is  constant  as 
warranted  by  Equation  6-22 . 

= constant;  on  cylinder  wall 

On  the  surface  of  the  sphere  the  stream  function  has  the 
same  value  as  the  centerline. 


'I^sphere  ® 
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Also  at  the  solid  surface  a no  slip  condition  exists  which 
makes  u = v = 0 

The  vorticity  is  calculated  by  using  its  definition.  A 
detailed  discussion  of  utilizing  different  differencing 
schemes  is  provided  by  Roache  (1982).  Regardless  of  the  wall 
orientation  or  the  boundary  value  of  i|r,  we  find 


rAn2 


0(  An) 


(6-23) 


where  An  is  the  distance  from  the  wall  point  (w)  to  the  next 
closest  vertical  point  (w+1).  With  this  definition  (Equation 
6-23)  there  is  no  problem  in  calculating  the  vorticity  for 
walls  aligned  with  the  cylindrical  grid.  However,  for  an 
inclined  surface,  such  as  that  presented  by  the  sphere 
surface,  some  interpolation  needs  to  be  done. 

For  the  primitive  variable  approach  a no-slip  condition 
results  in  u=v=0  for  the  solid  walls. 

Outflow.  The  outflow  is  labelled  as  surface  (4).  It  is 
a challenging  problem  in  computational  fluid  dynamics  to 
choose  this  condition  beforehand.  It  is  highly  problem 
dependant.  For  this  problem  the  boundary  condition  for 
stream  function  and  vorticity  as  used  by  Fanning  and  Mueller 
(1973)  was  utilized.  This  means  that  a linear  extrapolation 
of  values  in  the  axial  direction  was  made. 

“ 2 “ %M-2,j  (6-24) 

>IM,j  ~ 2 ” ClM-2,j 


(6-25) 
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In  the  above  equation  IM  is  the  x value  of  the  outflow 
boundary ‘This  condition  allows  the  flow  to  continue  in  the 
general  direction  as  dictated  by  the  two  interior  points 
adjacent  to  the  boundary.  The  velocities  at  the  outflow  were 
computed  using  centered  and  one-sided  differencing 


2rAr 


(6-26) 


IM,j 


rAz 


(6-27) 


By  using  the  above  method  stable  results  were  obtained.  It 
must  be  noted  that  'wiggles'  are  produced  if  proper  outflow 
conditions  are  not  used.  This  has  been  discussed  by  Roache 
(1982)  pg.  161.  Also  upwind  differencing  scheme  and  the  use 
of  gradient  condition  d^/dz  = 0 do  not  produce  wiggles  or 
oscillations  in  the  solutions. 

For  the  primitive  variable  approach,  linear  extrapolation  of 
velocity  was  used. 

Initial  condition.  Two  types  of  initial  conditions  were 
used.  The  first  was  the  extension  of  the  inflow  condition 
throughout  the  interior  region.  However,  for  the  stream 
function  vorticity  method  only  the  flat  velocity 
distribution  as  the  inflow  condition  converged  easily.  The 
second  method  utilized  a converged  solution  at  some 
different  Reynolds  number. 
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6.6  Calculation  of  Flow  Variables 
The  two  different  methods  used  to  calculate  the  flow 
variables  warrant  different  implementation  strategy. 

Streamf unction-Vorticity . The  following  steps  are 
undertaken  to  calculate  the  flow  variables  by  this  method. 

a)  The  grid  points  are  first  generated.  This  involves  the 
numerical  grid  generation  using  elliptic  equations  like 
Equations  6-13  and  6-14. 

b)  Once  the  grid  points  within  the  physical  space  are 
determined  an  initial  value  for  all  the  flow  variables  is 
assigned  to  them.  The  process  starts  by  providing  values  of 
the  streamfunction  and  correspondingly  computing  the  values 
of  velocity  and  vorticity  through  proper  relational 
equations. 

c)  The  vorticity  transport  equation  (Equation  6-7)  is  then 
solved  to  get  the  value  of  vorticity  at  the  next  time  step. 
This  now  gives  us  the  new  values  of  vorticity  at  each  grid 
point . 

d)  From  the  Poisson's  equation  (Equation  6-9),  the  value  of 
the  streamfunction  is  updated.  A point  successive  relaxation 
technique  was  used  to  solve  the  Poisson's  equation. 

e)  The  new  values  of  the  streamfunction  now  help  in  deriving 
the  velocity  components  in  the  'z'  and  'r'  directions 
respectively,  through  Equation  6-8. 
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f)  The  new  boundary  values  of  the  streamf unction  (Y)  and 
vorticity  ( $ ) are  calculated  by  using  new  i|;  and  5 values  at 
the  interior  points. 

g)  The  solution  proceeds  by  going  back  to  c)  in  calculating 
§ and  follows  steps  d)  through  f).  This  procedure  goes  on 
until  a convergence  level  or  the  desired  time  level  is 
reached  and  this  provides  the  final  solution. 

Primitive  Variable  Approach.  The  values  of  the  flow 
variables  for  this  method  are  obtained  by  using  the  SIMPLE 
algorithm.  This  algorithm  is  discussed  in  detail  by  Patankar 
(1980)  and  is  briefly  outlined  here. 

a)  As  in  the  case  of  the  streamfunction-vorticity  method, 
initially  the  grid  points  within  the  physical  space  are 
first  laid  out.  In  this  case  boundary  fitted  co-ordinate  are 
not  used  and  an  (r-z)  axisymmetric  cylindrical  system  is 
used. 

b)  The  initial  condition  is  given  by  providing  a value  of 
pressure  (p*)  at  each  grid  point. 

c)  The  momentxim  Equations  6-1  and  6-2  are  discretized  using 
the  control  volume  approach  and  following  Appendix  B 

^nb^nb'*'^'*'^e(PD~Pe)  28) 

^nb^nb*b-^An(Pp-Pn)  ( ) 

The  summation  sign  is  used  to  add  the  contribution  of  the 
four  (E,W,N,S)  neighboring  points.  'A'  is  the  cross 
sectional  area  and  'b'  is  the  source  term.  The  velocity 
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components  are  used  in  a staggered  form.  From  the  guessed 
value  of  the  pressure  field  (p  ) and  the  neighboring 
velocities  Uj^j^  and  a tentative  set  of  velocities  u*  and 

V*  can  be  estimated  from  Equations  6-28  and  6-29 


^nb^nb ( Pp  “Pe  ) /^e 
=E  ^nbKb  ( Pp  ~Pn  ) /^n 


d)  The  velocities  u*  and  v*  will  in  general  not  satisfy  the 
continuity  equation  (Equation6-3 ) . If  a correction  p'  is 
applied  to  the  guessed  value  of  p*  such  that 

p = p*  + p'  (6-32) 

then  one  gets  the  corresponding  corrected  value  of  the 
velocity  u'  which  is  defined  by 

u = u*  + u'  (6-33) 

e)  By  subtracting  Equation  6-30  from  Equation  6-28  one  gets 


E ^nb^nb’^^eiPp  Pe  ) 


(6-34) 


The  term  Lajj^^nb  considered  as  negligible  and  is 

dropped.  Equations  6-28,  6-30,6-32  and  6-33  combine  to  give 
a velocity  correction  formula 

Ug  = u*  + de(p^-p|.)  (6-35) 

where  dg  = Ag/ag.  In  a similar  manner  we  get  for  Vj^, 

Vn  = < + <^n(Pp-Ptf) 


where  d^  = A^/a^ 


(6-36) 
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f)  The  pressure  correction  term  (p')  is  derived  from  the 
continuity  equation  and  is  given  by  the  algebraic  form: 

^pPp  ~ ^ePe  ^wPw  ^nPn  ^sPs  ^ (6“37) 

where  a^  = (pft.d)g 

= (p^d)w 

a„  = (p^d)n 
ag  = (fAd)g 

b = (fAu*)„-(fAu*)3+(pft.u*)g-(jAu*)n-(f:t>  - pg)(Av/At) 
Here  ' p'  is  the  density  of  the  fluid  , 'AV'  is  the  volume  of 
the  control  volume  used  for  discretization  and  At  is  the 
increment  in  time. 

g)  The  new  pressure  term  p'  is  used  as  the  new  guessed 
pressure  p*.  Calculation  is  returned  to  step  'c'  and  the 
procedure  is  continued  until  the  desired  convergence  is 
reached. 


6 . 7 Test  Cases 

Before  starting  the  solution  to  our  problem,  some  test 
cases  were  run  to  determine  the  efficacy  of  the  coded 
programs.  The  deviation  of  the  computed  results  from  known 
results,  serves  to  analyze  the  effectiveness  of  the  computer 
code.  Different  cases  were  run  for  the  stresunfunction- 
vorticity  method  and  for  the  primitive  variable  method. 

In  the  streamfunction-vorticity  case  a boundary  fitted 
coordinate  system  was  used.  As  a preliminary  test,  fluid 
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flow  through  a cylindrical  pipe  was  considered.  It  is  well 
known  that  for  Icuminar  flow  a parabolic  profile  is  attained 
after  the  flow  has  fully  developed.  To  test  this,  a cylinder 
with  sufficient  length  was  considered  and  a flat  velocity 
profile  at  the  tube  inlet  was  taken.  The  flow  Reynolds 
nximber  based  on  this  mean  velocity  was  50.  The  grid  was 
generated  to  get  its  concentration  towards  the  centerline. 
The  transformation  of  the  physical  coordinate  (r-z)  to  the 
computational  coordinate  ($-T|)  is  given  by 

(6-38) 

(fl+l)-(Q-l){[(fl+l)/(Q-l)]^'^}  (6-39) 

[ (fi+l)/(Q-l) 

where  Q is  a clustering  parameter  within  the  range  of  1 to 
o.  As  the  value  of  approaches  1,  more  grid  points  are 
clustered  near  the  bottom  surface  and  which  happens  to  be 
the  centerline  of  the  cylinder  in  our  case.  Figure  6-8  shows 
the  grid  lines  as  they  appear  on  the  physical  and 
computational  space.  It  was  observed  that  the  value  of  the 
parameter  Q determined  the  accuracy  of  results  obtained. 

When  Q had  a value  of  2 the  error  in  determining  the 
centerline  velocity  was  0.35%.  With  the  decrease  in  the 
value  of  Q,  the  error  increases.  Figure  6-9  shows  the  flow 
parameters  as  obtained  by  using  P equal  to  1.6.  Figure  6-10 
is  a plot  of  the  absolute  error  obtained  in  calculating  the 
centerline  velocity  by  changing  the  nvimerical  value  of  Q. 
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a.  Physical  domain  (£2=1.02)  b.  Physical  domain  (Q=1.2) 


c . Computational  domain 


Figure  6-8  Physical  and  computational  domain  generated  by 

transformation  Equations  6-38  and  6-39. 


a.  Grid  lines 


b.  Stream  function  contour  lines 


c.  Vorticity  contour  lines 


d.  Velocity  vectors 


Figure  6-9.  Grid  lines  and  flow  parameters  for  (3=1.2 
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Figure  6-10.  Absolute  error  in  determining  the  centerline 

velocity  for  Reynolds  number  of  50 
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As  a second  case,  a cylindrical  pipe  was  again 
considered,  but  it  had  a bell  shaped  constriction.  Figure  6- 
11  a.  shows  the  shape  of  the  physical  space.  Lee  and  Fung 
(1970)  had  perforaed  some  numerical  study  on  such  a system 
and  Bentz  and  Evans  (1975)  conducted  an  experimental  study 
by  using  a laser  Doppler  anemometer.  The  axisymmetric 
constriction  was  specified  by  the  equation 

Dn  1 2 ^ 

r=  ° [ l-iexp(  -16£^  ) ] ( 6-40 ) 

were  'r'  is  the  value  of  the  local  radius  of  curvature  of 
the  wall,  Dq  is  the  diameter  of  the  tube,  and  'z'  is  the 
distance  in  the  axial  direction.  Reynolds  numbers  of  2 and 
25  were  used  to  compare  results  with  the  published  values. 
The  results  are  plotted  and  shown  in  Figure  6-11  b and  6-11 
c.  It  can  be  noticed  that  the  error  from  the  present  work  is 
larger  than  the  previous  numerical  work.  Lee  and  Fung  (1970) 
had  however  used  confomal  mapping  to  generate  the  grid 
lines  and  this  eliminated  any  error  arising  during  the 
mapping  from  physical  to  computational  regime.  Also,  it  was 
observed  in  the  present  work  that  the  error  propagated 
during  the  iterative  computation  of  the  streamfunction 
values  and  not  during  the  calculation  of  the  vorticity 
transport  equation.  This  was  tested  by  using  analytical 
values  of  metrices  during  the  iterative  calculation  of  the 
streamfunction.  The  vorticity  values  obtained  from  the 
transport  equation  did  not  vary  much  if  numerically 
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a.  Axisymmetric  bell  shaped  constriction 


b.  Reynolds  number  = 2 


c . Reynolds  number  = 25 


Figure  6-1 1 Grid  points  and  nondimensional  velocity  values 

at  the  center  plane 
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generated  values  of  the  metrices  were  used  instead  of  the 
analytically  known  values.  However  this  point  can  be  stated 
emphatically  only  after  rigorous  testing  of  several  cases 
and  which  was  not  done  here. 

When  the  streamfunction-vorticity  method  was  applied  to 
problem  of  sphere  within  a cylinder,  convergence  became  a 
major  problem.  This  can  be  explained  in  the  light  of  the 
test  cases  run.  In  order  to  get  good  values  near  the  surface 
of  the  sphere,  the  value  of  the  stretching  factor  P is 
chosen  close  to  1.  But  as  can  be  seen  from  Figure  6-10,  even 
in  simple  cases  this  resulted  in  significant  error.  Also  as 
had  been  explained  earlier,  when  the  flow  is  not  aligned  to 
the  grid,  error  arising  from  numerical  diffusion  becomes 
dominant.  In  the  present  geometry  both  these  factors  play  a 
contributory  role.  In  the  case  of  bell  shaped  constriction, 
the  value  of  the  Reynolds  number  was  less  and  convergence 
did  not  create  a problem  although  the  results  were  not  quite 
accurate.  Similarly, in  the  case  of  a sphere  within  a 
cylinder,  the  convergence  was  not  difficult  when  flow 
Reynolds  number  (based  on  the  tube  dicuaeter  and  mean  flow 
velocity)  was  equal  to  30  and  55  respectively.  If  the 
Reynolds  number  is  based  on  the  sphere  diameter  then  in  both 
the  cases  it  did  not  exceed  20.  Figure  6-12  shows  some  of 
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a . p=  0.2 


b.  p=0.35 


c.  P=0.5 


Figure  6-12  Different  diameter  ratio  (d/D)  considered 
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the  different  sphere  to  cylinder  diameter  ratios  ( P)  that 
were  considered.  Figures  6-13  and  6-14  give  an  estimate  of 
how  the  streamf unction,  vorticity  and  the  velocity  vectors 
look  like  for  different  Reynolds  numbers.  Also  from  Figure 
6-15,  one  finds  that  the  value  of  vorticity  on  the  surface 
of  the  sphere  did  not  fall  below  zero.  This  indicates  that 
the  flow  did  not  separate  for  the  Reynolds  number 
considered.  This  is  fine  since  it  has  been  found  that  the 
flow  separation  takes  place  when  the  Reynolds  number  (based 
on  sphere  diameter)  reaches  20.  Beyond  Reynolds  number  of 
100  convergence  was  not  obtained.  Some  of  the  possible 
reasons  for  this  have  been  discussed  earlier. 

As  a result  numerical  testing  was  carried  out  by 
primitive  variable  approach.  This  program  incorporated  a 
linear  flux  spline  scheme  to  reduce  error  when  the  flow  is 
not  aligned  with  the  grid  lines.  Although  the  surface  of  the 
sphere  was  not  traced  out  correctly  through  a 'r-z' 
coordinate,  but  results  obtained  were  meaningful.  The 
robustness  and  exactness  of  the  program  had  already  been 
tested  by  Karki  et  al.  (1988),  and  no  further  test  cases  are 
shown  here.  Figure  6-16  shows  the  separation  angle  as  a 
function  of  the  sphere  to  cylinder  diameter  ratio.  The  flow 
Reynolds  number  was  taken  as  100  and  the  result  shows  an 
interesting  similarity  with  the  results  obtained  through 
boundary  layer  analysis.  The  separation  angle  reaches  a 
minimum  when  the  diameter  ratio  is  around  0.5.  Since 
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different  flow  rates  are  used  while  comparing  the  numerical 
results  with  the  boundary  layer  analysis,  the  magnitude  of 
the  separation  angle  differs.  In  boundary  layer  analysis  the 
Reynolds  number  was  around  3000,  but  in  the  numerical  case 
it  is  much  less  (to  consider  an  axisymmetric  flow).  Figure 
6-17  gives  a plot  of  the  separation  angle  as  a function  of 
Reynolds  number  with  a constant  sphere  and  cylinder  size. 
This  has  also  been  compared  with  the  flow  over  an  unbounded 
sphere.  In  the  case  of  free  flow,  the  separation  angle  is 
given  by  0g  = 180  - 42.5[ln(Re/20) but  in  the  bounded 
case  no  explicit  equation  has  been  proposed  here  due  to  the 
lack  of  several  data  points.  The  drag  force  as  obtained  from 
the  experiment  has  been  calculated  with  the  numerical 
program.  Since  the  sphere  surface  was  not  properly  obtained 
from  the  'r-z'  coordinates,  a good  result  should  not  be 
anticipated.  Figure  6-18  provides  a comparison  between  the 
experimental  and  numerical  results.  By  increasing  the  number 
of  grid  points  on  the  sphere  surface,  the  drag  force  results 
initially  approach  the  experimental  values.  But  further 
increase  in  the  grid  point  nmnber  did  not  help  the  cause  to 
a large  extent.  Finally  some  prdicted  flow  patterns  for 
P( ratio  of  diameter  between  sphere  and  cylinder)  of  0.5  and 
0.1  are  shown  in  Figures  6-19  through  6-22  for  Reynolds 
numbers  of  50  and  250. 
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b.  Vorticity  (Re=30) 


c.  Velocity  vector  (Re=30) 

Figure  6-13.  Flow  patters  for  Reynolds  number  of  30 


a.  Streamline  (Re=55) 


b.  Vorticity  (Re=55) 


c.  Velocity  vector  (Re=55) 

Figure  6-14  Flow  patters  for  Reynolds  number  of  55 
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(3=0.35) 
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Plot  of  separation  angle  and  different  sphere 
to  cylinder  diameter  ratio 


Figure  6-16 
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Figure  6-17  Plot  of  separation  angle  for  different  flow 
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Figure  6-18  Drag  force  calculation  on  sphere  (6=0. 
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a.  Reynolds  nuniber  = 50,  ^,epar«tion=  165.08  degrees 


b.  Reynolds  number  =»  250,  ^separation®  122.86  degrees 
Figure  6-19  Velocity  vectors  for  P=0.5 
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a.  Reynolds  number  = 50,  ^separation=  165.08  degrees 


b.  Reynolds  number  = 250,  <I>.eparation=  122,86  degrees 
Figure  6-20  Velocity  (U)  profile  for  |3=0.5 
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a.  Reynolds  nxunber  = 50,  4>,eparation=  165.08  degrees 


b.  Reynolds  number  = 250,  ^separation=  122.86  degrees 
Figure  6-21  Pressure  contours  for  |3=0.5 
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Reynolds  number  =50  b.  Reynolds  number  =250 

^separation*  173.95  degrees  ^separation*  142.76  degrees 

Figure  6-22  Velocity  profile  for  p=0.1 
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a.  Reynolds  number  =50  b.  Reynolds  number  =250 

0.epar.tion=  173.95  degrees  Q.ep*rauon=  142.76  degrees 


Figure  6-22  Pressure  contours  for  (3=0.1 
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CHAPTER  7 

SUMMARY  AND  CONCLUSION 

The  current  study  tries  to  further  the  investigation  of 
a fundamental  fluid  flow  problem;  the  drag  force  measurement 
on  a cylindrically  bounded  sphere.  This  work  is  then 
extended  for  the  viscosity  measurements  of  gases.  The 
methodology  proposed  for  such  viscometric  measurements  has 
very  attractive  features  when  the  gas  is  corrosive  or  has  a 
high  temperature. 

Previous  work  in  the  area  of  drag  measurement  for 
spheres  in  a cylinder  had  focussed  on  sedimenting  spheres. 
For  very  low  flow  rates  analytical  work  was  performed  by 
neglecting  the  convective  terms.  Later  improvements  were 
initiated  by  taking  other  factors  into  consideration-  like 
linearizing  the  inertial  term,  taking  the  wall  effect  etc. 
The  results  of  these  studies  provided  a means  of  measuring 
the  viscosity  of  fluids.  However,  since  gases  have  a low 
viscosity  the  sedimentation  study  had  a practical 
implication  for  liquid  viscosity  measurements  only.  Among 
other  ideas,  the  capillary  tube  and  the  oscillating  body 
were  of  less  restrictive  nature  and  as  such  are  used  as  a 
primary  viscometer  till  date.  However  even  these  primary 
instruments  run  into  difficulty  when  the  fluid  is  of 
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corrosive  nature.  Unresearched  materials  compatibility  and 
the  unknown  performance  of  the  measuring  device  give  rise  to 
the  uncertainty  in  the  final  results.  The  present  work 
involves  the  drag  measurement  on  a sphere  which  is  rigidly 
supported  at  the  center  of  a cylinder.  The  viscosity  is 
ascertained  by  knowing  an  accurate  plot  of  the  coefficient 
of  drag  and  the  flow  Reynolds. number . For  an  unknown  gas  the 
coefficient  of  drag  can  be  had  from  its  flow  rate  and  then 
from  the  plot  the  Reynolds  number  is  obtained.  This  provides 
the  numerical  value  of  the  viscosity.  An  interesting  feature 
of  the  setup  is  that  the  measuring  devise  is  not  exposed  to 
the  fluid.  As  such  during  the  measurements,  one  is  not 
concerned  with  the  reaction  between  the  fluid  and  the  drag 

I 

measuring  device.  This  helps  the  methodology  to  be 
successfully  applied  for  any  corrosive  gas.  For  high 
temperature,  the  sphere  is  only  heated,  and  the  temeperature 
gradient  is  felt  near  the  sphere  surface  only.  As  a result 
the  flowing  gas  has  a much  lower  temperature  as  it  exits  the 
system  and  most  the  problems  related  with  the  high 
temperature  is  eliminated.  In  the  actual  study,  only  clean 
gases  like  Oxygen  and  Helium  were  used  at  room  temperature, 
and  the  results  were  quite  good. 

Very  low  flow  rates  are  however  desirable  in  order  to 
determine  an  accurate  value  of  viscosity.  At  creeping  zone, 
Haberman  and  Sayre  (1958)  give  the  drag  force  on  a fixed  and 
rigid  sphere  in  a Poiseuille  flow  as 
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DRAG  = 6 n^R  U 


max 


K 


(7-1) 


where 


K= 


1-2.3^-0.202173^ 

3 


1-2. 105/3+2.0865)3^-1.7068/3^+0.72603/3^ 


J3  is  the  diameter  ratio  between  the  sphere  and  the  cylinder, 
jU  is  the  viscosity  of  the  fluid,  R is  the  radius  of  the 
sphere  and  is  the  centerline  velocity  for  Poiseuille 

flow.  In  order  to  get  the  maximum  sensitivity  in  the  final 
results,  one  needs  to  choose  a large  sphere.  For  the  present 
microbalance  the  maximum  weight  restriction  is  3.5  grams.  If 
one  uses  teflon  ball  then  we  can  choose  a sphere  of  diameter 
around  12.5  mm.  and  for  nitrogen  the  drag  forces  are: 


Ratio  (J3) 

Draa  Force  /ka) 
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3.3484 

E-6 
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• 
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7.0558 

E-6 

o 

• 

VO 

1.5959 

E-5 

Thus  the  drag  force  can  be  measured  with  a resolution  of  one 
part  per  thousand  (since  the  sensitivity  of  instrument  is 
0.1  fug  and  the  drag  force  is  in  mg).  Instrximent  this  precise 
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were  not  available  earlier,  and  gas  viscosity  measurements 
by  direct  drag  measurement  were  previously  ignored. 

The  sphere  is  suspended  by  a nylon  string  and  it  can 
have  its  presence  felt  in  two  ways.  Firstly  the  drag 
measured  in  the  experiment  is  the  sum  of  the  drag  on  the 
sphere  and  the  string.  In  order  to  eliminate  the  drag  on  the 
string,  two  different  string  lengths  for  the  same  geometry 
and  flow  conditions  were  used.  This  method  provided  very 
accurate  results  as  is  seen  from  Table  3-1.  Secondly  the 
fluid  flow  parameters  for  an  isolated  sphere  versus  a sphere 
supported  by  a string  can  be  different.  However  it  has  been 
found  by  several  authors  that  if  the  diameter  of  the  string 
is  small  in  comparison  to  that  of  the  sphere  then  the  flow 
remains  unaffected.  Seeley  (1973)  found  the  flow  pattern 
affected  when  D/d  was  equal  to  18.8  and  10.9,  where  D was 
the  sphere  diameter  and  d was  the  string  diameter.  Taneda 
(1956)  found  that  if  D/d  ratio  is  30  then  the  wake  remained 
unaffected  from  Reynolds  number  of  5 to  300.  Also  Kalra  and 
Uhlher  (1971)  investigating  the  flow  over  a supported  sphere 
concluded  that  the  string  had  negligible  effect  on  the 
separation  angle  and  the  wake  dimension  when  D/d  was  36  and 
the  Reynolds  number  ranged  from  30  to  750.  In  our  experiment 
the  string  was  0.07  mm  in  diameter  and  the  sphere  diameter 
ranged  from  7.88mm  to  12.67mm  and  as  such  its  effect  on  the 
flow  can  be  safely  neglected. 
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One  of  the  parameters  of  paramount  importance  for  the 
fluid  flow  over  a bluff  body  is  the  separation  angle.  It  is 
very  difficult  to  ascertain  it  correctly  by  experimentation 

j 

or  by  numerical  approach.  In  this  study  a boundary  layer 
analysis  is  performed  to  observe  the  dependence  of  the 
separation  angle  on  the  diameter  ratio  between  the  sphere 
and  the  cylinder  (B).  This  study  is  based  on  the  assumption 
that  the  Reynolds  number  is  fairly  high  but  low  enough  to 
have  a laminar  boundary  layer.  The  results  show  that  the 
separation  angle  reaches  a minimum  around  J3  = 0.5. 
Interestingly  enough  the  same  pattern  is  also  provided  by 
the  numerical  studies.  However  the  Reynolds  number  for  the 
numerical  studies  should  be  low  so  that  the  flow  is  two 
dimensional  throughout.  This  is  necessary  since  the  code  was 
written  for  axisymmetric  case  only. 

The  results  of  the  drag  force  calculation  is  best 
modelled  by  obtaining  a relationship  between  the  coefficient 
of  drag  (Cq)  and  the  flow  Reynolds  number.  Previous  authors 
who  had  performed  investigation  had  come  up  with 
correlations  linking  the  coefficient  of  drag  for  an 
unbounded  sphere  (C^q,)  to  the  coefficient  of  drag  for  a 
bounded  sphere  (C^)  for  the  scune  Reynolds  number.  However 
different  correlations  were  propounded  due  to  different 
results  obtained  in  the  experiments.  One  of  the  disturbing 
features  has  been  the  nebulous  definition  of  the  Reynolds 
number  while  comparing  the  unbounded  and  the  bounded  sphere. 
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Some  authors  have  used  the  mean  velocity  in  the  tube  and  the 
others  have  chosen  the  centerline  velocity  as  the 
representative  velocity  to  define  the  Reynolds  number.  There 
has  also  been  a dilemma  over  the  choice  of  cylinder  diameter 
over  the  sphere  diameter  in  this  regard.  As  far  as  the 
velocity  is  concerned,  a problem  arises  since  there  is  a 
velocity  distribution  for  a bounded  sphere  but  an  unbounded 
sphere  is  assvimed  to  encounter  a constant  mean  velocity. 

Thus  based  on  the  physics  of  the  problem  the  current  work 
proposes  a characteristic  velocity  which  gives  the  same 
dynamic  force  as  exerted  by  the  velocity  distribution  seen 
by  the  sphere.  It  is  given  by  V {l-(2J3^/3)  + (J3^/5) } 

where  is  the  centerline  velocity,  and  B is  the  diameter 

ratio  between  the  sphere  and  the  cylinder.  In  the  above 
relationship  the  mean  is  obtained  in  a linear  way.  If  an 
area  averaging  is  done  then  the  characteristic  velocity  is 
given  by  Reynolds 

number  is  based  on  these  characteristic  velocities  and  the 
diameter  of  the  sphere,  then  a simple  relationship  is 
obtained,  as  given  by  Equation  5-18  or  5-19,  for  the 
coefficient  of  drag  between  a bounded  and  unbounded  flow. 

Finally  an  attempt  was  made  to  explain  the  fluid  flow 
phenomenon  by  numerical  studies.  The  two  popular  approaches 
to  incompressible  fluid  flow  are  the  streamfunction- 
vorticity  method  and  the  primitive  variable  method.  During 
streamfunction-vorticity  method  a boundary  fitted  coordinate 
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system  was  applied.  Although  it  appears  to  be  attractive  for 
the  geometry  in  hand,  but  the  results  can  be  erroneous.  The 
grid  was  generated  by  elliptic  grid  generation  technique. 
This  resulted  in  smooth  grid  lines,  but  relatively  sparse 
grid  density  close  to  the  sphere  surface.  As  has  been  found 
in  test  cases,  for  boundary  fitted  coordinate  in 
streamfunction-vorticity  method,  the  convergence  did  not 
guarantee  a good  result.  A stretching  parameter  of  1.2  was 
used  in  the  present  case.  The  root  cause  of  the  deviation  in 
the  answers,  from  exact  values,  had  been  during  the 
iterative  solution  for  the  streamf unction.  This  has  been 
discussed  in  chapter  6 and  concluded  by  using  analytical 
values  for  the  metrices.  The  program  also  showed  little 
signs  of  convergence  after  the  Reynolds  number  reached  100. 
Attention  was  paid  to  the  primitive  variable  method  which 
employed  r-z  coordinates.  Here  the  sphere  surface  was  not 
properly  described  by  the  grid  lines,  but  the  results  were 
satisfactory.  In  this  latter  method,  a linear  flux  spline 
was  used  to  model  the  fluxes  in  the  control  volume,  during 
the  process  of  discretization.  This  has  been  found  to  handle 
problems  in  an  efficient  way  when  the  flow  is  highly  non- 
aligned  to  the  grid  lines.  The  results  obtained  show  the 
same  trend  as  has  been  shown  by  the  experiment.  From  the 
experience  gathered  in  this  process,  the  following 
recommendations  can  be  made  for  future  studies.  Firstly  a 
body  fitted  system  for  the  primitive  variable  approach 
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should  be  used.  By  using  different  stretching  factors,  the 
results  will  not  change  as  dramatically  as  in  the 
streamfunction-vorticity  method.  As  had  been  found  in  the 
streamfunction-vorticity  method,  during  the  iteration  for 
the  vorticity  and  streamf unction,  the  niimerical  results  for 
the  metrices  gave  a satisfactory  value  for  the  vorticity  but 
not  for  the  streamf unction.  This  situation  does  not  exist  in 
the  primitive  variable  method,  were  conservative  equations 
like  those  for  the  computation  of  vorticity  is  used  all  the 
time.  However  one  will  still  not  get  good  resolution  for  the 
front  and  aft  of  the  sphere,  by  using  elliptic  equations.  If 
one  uses  a grid  that  is  forced  to  follow  the  sphere  surface, 
then  the  value  of  the  metrices  may  cause  the  results  to 
diverge  from  the  actual  results.  One  of  the  remedies  to 
overcome  this  is  to  use  a multiblock  composite 
configuration.  In  this  procedure,  the  physical  space  is 
divided  up  into  convenient  number  of  regions  (blocks)  and  a 
body  fitted  coordinate  is  used  for  each  of  them.  Caution 
must  be  taken  to  link  each  of  the  blocks  through  the  common 
layer  they  include  amongst  themselves.  Figure  7-1  shows  the 
mapping  of  the  physical  space  to  the  computational  space  by 
this  procedure. 
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b.  Computational  region 


Figure  7-1  Multiblock  composite  configuration 
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NUMERICAL  VALUES  USED  FOR  BOUNDARY  LAYER  CALCULATION 

The  calculations  were  carried  out  by  expanding  the 
outer  velocity  (U),  the  local  radius  of  curvature  (r)  and 
the  stream  function  (7)  respectively  as  a power  series 
(Schlichting  1968): 


u = 

U^x 

uax^ 
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(A- 
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riX 
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where  x is  measured  along  the  body  as  shown  in  Figure  4-2. 
The  terms  f^^,  ±2/  ^5/  ®tc.  as  needed  in  Section  4.4  have  to 
be  calculated.  For  this,  two  sets  of  coefficients  need  to  be 
ascertained.  The  first  is  u^,  U3,  Ug  ...  etc.  whose  values 
have  already  been  obtained  in  Section  4.4  and  are  repeated 
here  for  convenience: 

Ui=  K UgTj^/r  , U3=  K UoT3/r^  , Ug=  K UoTg/r®'  . . . and  so  on. 
Here  K is  a constant,  Uq  is  the  velocity  hitting  the  sphere 
at  the  front  stagnation  point,  Tj^,T3,T5, . . . , etc.  are  the 
coefficients  obtained  by  expanding  the  term  (sin$)/(l- 
6^sin^$)  and  r is  the  local  radius  of  curvature  of  the 
sphere.  The  second  set  of  coefficients  is  r3^,r3,rg,  ...  etc. 
For  a sphere  the  local  radius  of  curvature  (r)  is  given  by 
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Rsin$,  where  $ is  the  angle  subtended  by  the  length  x of  the 
sphere  of  radius  R at  the  center. 

Therefore,  r = Rsin$ 

= R($  - $^/3!+  $^/5!  - $'^/7I  + ...  ) 

= R(X/R  - X^/(R^3!)  + X^/(R®5!)  - x"^/ ( r'^7  ! ) + . . . 

(A-4) 

Comparing  coefficients  of  Equations  (A-2)  and  (A-4)  we  get: 
r^=  1,  r3=  1/(6R^),  rg=  1/(120R^),  Vj=  1/(5040r'^)  ... 

The  accuracy  of  this  power  series  type  of  solution  depends 
on  the  number  of  terms  taken.  Here  the  expansion  was  carried 
out  for  seven  terms  and  following  Scholkemeier  (1949)  we 
get: 
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i.e. 
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In  order  to  calcualte  the  shear  stress  the  values  at  the 
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These  can  be  sustituted  in  the  following  equations: 
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1 1 
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(A-10) 


to  get  the  needed  results. 


APPENDIX  B 


DISCRETIZATION  USING  CONTROL  VOLUME 


The  generalized  time  dependant  governing  equation  in 
axisymmetric  cylindrical  coordinate  (r-z  with  symmetry  about 
0)  that  needs  to  be  discretized  so  that  it  reduces  to  an 
algebraic  form  is 

d(p<p)  4.  ^ < P^r*P ) 4.  ^ ^ d d ,j,d(p  d , T(p 

—Wt JE — 


(B-1) 


where  t is  the  time,  V^.  and  are  the  velocities  in  the  r 
and  z direction  respectively,  p is  the  density  of  the  fluid, 
$ is  a dependent  variable  and  r is  the  diffusion 
coefficient. 

Figure  B-1  shows  a typical  control  volume  whose  space 
satisfies  Equation  B-1.  Integrating  each  term  of  the 
governing  equation  with  respect  to  the  dz,dr  and  dt  we  get 
for  the  time  dependent  term  (first  term  of  Equation  B-1); 


{t+At)ne 

i n 


g(p0) 


dzdrdt  = { Pp(pp-pp<t>p)  AzAr 


(B-2) 


where  the  index  '0'  denotes  value  at  time  t. 

The  second  and  third  term  on  the  left  side  of  Equation 
B-1  are  the  convective  terms  and  their  integration  yields; 
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(t+At)ne 

I n 


d ( pV-d) ) 

dzdrdt  = [ ( pVr<l> ) 


n 


( pVr<p ) ^ ] ^zAt 


(B-3) 


(PVr)n  <P^r) 


}*Pp 


ipVr) 


(pg  ] Az  At 


In  Equation  B-3  the  value  of  $ at  the  n and  s face  is 
obtained  by  averaging  i.e. 

$n=($P+«N)/2  and  $g=($p+$s)/2 
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Figure  B-1  A Typical  Control  Volume  for  Discretization 


Similarly  for  the  other  convective  term  we  get 
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(t+At)ne 

I II 


d ( pV^<l> ) 


— dzdrdt  = [ ( pV^(p ) ^ 
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w 


}*pp 


ipv^) 


-<p„]^At 


The  first  two  terms  on  the  right  side  of  Equation  B-1  are 
the  diffusive  terras  and  we  can  integrate  them  as 


( t +At  )ne 


d<p 


d(f> 


_-r 


„ )ArAt 


(^)w 
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(B-5) 
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= rp  <^P_p 


= [ — ■■.g..-(^y+  . ^ (ps~{  - ^.9. + . ^ }(f)p]AzAt 

{Ax)n^^  {Ax)s  ^ (Ax)„  (Ax)s 

Finally  one  can  integrate  the  source  terra,  which  is  the  last 
term  of  Equation  B-1,  and  get 


{t*At)ne 

i n 


^{H)dzdrdt  = [(^)  -{^)  ]ArAt 
dr  r r „ r » 


(B-7) 
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We  can  observe  that  the  integrated  value  of  the  source  term 


can  be  expressed  as  S = S^,  + Sp$p.  Here  S^,  is  the  constant 
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term  independent  of  $p  and  Sp  is  the  coefficient  of  $p.  As 
has  been  discussed  by  Patankar  (1980)  for  convergence  by 
iterative  process  Sp  must  have  a negative  value.  If  the 
diffusion  coefficient  has  a constant  value  (i.e.  r=rjj=Pg) 
then  since  rg  is  always  less  than  r^^,  Sp  has  a negative 
value. 

Introducing  the  following  definitions: 


c = 

(pi)  As 

Ce  = 

(f»^z)e 

Ar 

Cw  = 

Ar 

Cn  = 

Az 

n 

CD 

III 

(P*^r)s 

Az 

we  can  substitute  these  terms 


; D = (r/Ax)  As 

; Dg  H (Pg/Axg)  Ar 

; D„  = (r„/Ax^)  Ar 

; D„  = (r„/Axn)  Az 

; Dg  = (Tg/AXg)  Az 

and  the  different  integral 


values  obtained  for  the  governing  Equation  B-1  to  get 


PpCpp  Pp<Pp 
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I.e. 


^f^£._^AzAr+(  ^ 
At  2 


W 
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+Dq+D„+D^+Ds)<Pp 


={d^-^)(Pe+{d„ 


(B-9) 

If  we  finally  provide  the  following  definition  as  given  by 
Patankar  (1980) 
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ag  = Dg  - Cg/2 

+ C„/2 

% “ “ ^n/2 

aj  = (pg  Az  Ar)/(At) 
b = Sj,  Az  + ap 

we  get  a = a^  + a^^  + ajj  + ag  + ap  - Sp  Az 


APPENDIX  C 


DISCRETIZATION  PROCEDURE  FOR  BOUNDARY  FITTED  COORDINATE 

In  the  boundary  fitted  coordinate  system  the  physical 
region  is  mapped  onto  a computational  space  which  has  a 
regular  geometry  with  evenly  spaced  grid  points.  Let  (5/11) 
be  the  variables  of  the  grid  points  in  the  computational 
domain  with  (r,z)  being  the  corresponding  variables  in  the 
physical  domain  as  described  by  the  axisymmetric  cylindrical 
system.  The  target  is  to  convert  the  generalized  governing 
Equation  like  Equation  6-10  or  Equation  B-1  from  (r,z) 
coordinates  to  ( 5 / t|)  coordinate  which  can  then  be  solved  in 
the  computational  region. 


Let  5 = 5(r,z) 

(C-1) 

T|  = T|(r,Z) 

(C-2) 

(C-3) 

II 

H 

+ 

(C-4) 

_p2 

3^2  ’’'3^2 

(C-5) 
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02  _j:2  02  „ 02  ^^2  02 

(C-6) 

Also  $2  ~ J 

(C-7) 

5r  — " Zt| 

(C-8) 

= - J rg 

(C-9) 

T^.  = J Ze 

(C-10) 

1 

where  J = 

(C-11) 

The  metrics  (C-7)  through  (C-11)  are  obtained  through  a 
procedure  found  in  most  computational  fluid  dynamics  text 
books.  The  Jacobian,  J,  is  interpreted  as  the  ratio  of  the 
areas  (for  two  dimensional  cases)  in  the  physical  space  to 
that  of  the  computational  space. 

With  these  values,  the  convective  terms  in  the 
governing  Equation  (B-1)  becomes 

z<P ) ^ ) 

^ (C-12) 

( P'^z'^^Tl  ) -JpVz4>  ( ^T1  ) 

( pv^<t>r^ ) +Jpv^ct> S(r^) 


d 

■Hi 


( pv^<f> ) =J 


d 


(ry^pv^<f>) 


(riPVz<P) 


(C-13) 


Similarly  ^ ( pV^<^)  (z^py^.^) 


(C-14) 


Let  T represent  time  during  calculations  in  the 
computational  space.  Noting  that  t = t one  has 
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d _ d c 3 d 

'Wt  ^<=‘9^ 

(C-15) 

But  = 0 and  = 1,  which  transforms  the  time 

dependant  terra  of  the  governing  equation  to  the 
computational  space  as 


(C-16) 

Taking  these  calculations  into  account  the  left  side  of 
the  governing  Equation  (B-1)  as  given  by 


-^  ( P^)  ( pv^<^)  ( 9VA) 

(C-17) 

is  now  transformed  for  the  computational  space  to 

^ { p0  ( r^V^-z^Vr) } { p^  ( z^Vr-r^V^ ) } 


If  we  introduce 

(C-18) 

U = r V - z V 
^ ■^T|''z  ^ryr 

(C-19) 

and  V = zsV^  - reV^ 
Equation  (C-18)  becomes 

(C-20) 

^ ( p^ ) ( m ) +J-4  ( ) 

(C-21) 

We  are  now  left  with  the  diffusive  and  source  terms.  The 
diffusive  terms  are  transformed  as  follows: 
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0(^  . _ 


=J  ^ (r 


3ri 


(C-22) 


(C-23) 


But,  J^(r,r||).J^(r,r(5,||*,i,||)> 

=J^(r,r(Jr,||-Jr5|^)) 

=J^{TJ{r^<f>^-r^r^<l>^)} 


(C-24) 


Similarly,  (j-{r|| ) -J  ® (rj{r{r^^{-r|^^) ) (C-25) 


Hence,  ^ (r|| ) -J^(rj(i-2^5-r;r^^,) } 


-J 


^{TJ(r^rr^<P^-rlcp^)} 


(C-26) 


Carrying  out  calculations  on  the  same  lines,  the  second 
diffusive  terms  becomes 


^ ‘^1?  > ) > 


(C-27) 


+J 


-^{JT(Z5>^-z^Z505)} 


Therefore  the  first  two  terms  on  the  right  side  of 
Equation  (B-1),  that  happen  to  carry  the  diffusive  flux,  are 


now  transformed  as 
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a 

~d^ 


[ rJ{  ( r| +z| ) (^^ - ( +2^Z5  ) } ] 

=J^  [rJ(Ai<^5-A2<^T,)  ] +J-4  [TJ(  A30^-A2<^5  ) ] 


where 

^2  “ ^T|^C  ^ ^T|^5 
A3  = z?  + r| 

Finally  the  source  term  is  transformed  as  follows: 


(C-28) 

(C-29) 

(C-30) 

(C-31) 


a 

■a? 


1 Try  ^ ^ ’*'*^^5  ^ ^ 


dz. 


=-,r  5 ) .,T  r<#> -T, 


(C-32) 


Therefore  by  taking  Equations  C-21,  C-28  and  C-32  we  get  the 
governing  equation  in  the  computational  plane  as 


a 

■a? 


( p<p ) +J-^  ( pu<p ) +J^  ( pv<p ) =J^  [ rj(  Ai^5  -A2<^t,  ) ] 

3 d ZeTd)  d ^n^<P 


(C-33 


Since  the  Jacobian,  J,  does  not  vary  with  time  t we  finally 
get 
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In  Equation  C-34  the  terms  U,V,A2,A2,A3  and  J are 
respectively  given  by  Equations  C-19,C-20,C-29,C-30,C-31  and 
C-11,  respectively. 

Discretization  of  Equation  (C-34)  is  performed  in  a 
similar  manner  as  shown  in  Appendix  A.  Integrating  over 
drd^dri  one  gets: 


0.0 


Pp<pp  Pp<pp  AEAti  ^3 

r. zl£  + (-^-^+_£--J+Pg  +p^  +p^  +p , 


)^p 


(C-35) 


In  Equation  (C-35)  the  various  terms  are  as  follows: 


C..  = 


Cw  = 


Cn  = 
Ca  = 


(fU)eATJ 

(fU)„AT1 


Dg  = (rJAi/A5)eATl 
D„  = (rJAi/A§)„Atl 
D„  = (rJA3/Ari)eA5 
Dn  = (rJA3/Ari)3A5 


and  the  source  term  S is  given  by: 


<f> 


d<p 


<P 


s=[(z^rx+rjA2-^)  ~(z^t:l+tja2:^)  ]Ati 


d<f> 
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* [{Z^T±+TJA2^)  ~(Z^T±+TJA2^)  ]A5 


(C-36) 
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Figure  C-1  A typical  control  volume  in  computational  space 


Equation  C-35  can  be  finally  cast  in  an  algebraic  form 


ap$p 


where  a™  = 


De  - Ce/2 
Dw  + C„/2 


(C-37) 


~ ^a/2 


ag  = ( ptA5Ari)/(JpAt) 


b = S + ap$p 
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It  is  interesting  to  note  that  during  calculation  the  point 
P is  influenced  not  only  by  points  E,W,N,and  S but  also  by 
points  NE,NW,SE  and  SW.  The  extra  four  points  NE,NW,SE  and 
SW  come  into  picture  while  calculating  the  source  term  S 
(Figure  C-1). 
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